source: trunk/ADOL-C/src/sparse/sparsedrivers.cpp @ 192

Last change on this file since 192 was 192, checked in by awalther, 9 years ago

check in patch for deep copy for set_HP and get_HP

  • Property svn:keywords set to Author Date Id Revision
File size: 36.8 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     sparse/sparsedrivers.cpp
4 Revision: $Id: sparsedrivers.cpp 192 2011-01-17 20:05:16Z awalther $
5 Contents: "Easy To Use" C++ interfaces of SPARSE package
6 
7 Copyright (c) Andrea Walther
8 
9 This file is part of ADOL-C. This software is provided as open source.
10 Any use, reproduction, or distribution of the software constitutes
11 recipient's acceptance of the terms of the accompanying license file. 
12 
13----------------------------------------------------------------------------*/
14#include <adolc/sparse/sparsedrivers.h>
15#include <adolc/oplate.h>
16#include <adolc/adalloc.h>
17#include <adolc/interfaces.h>
18#include "taping_p.h"
19
20#if defined(ADOLC_INTERNAL)
21#    if HAVE_CONFIG_H
22#        include "config.h"
23#    endif
24#endif
25
26#if HAVE_LIBCOLPACK
27#include "ColPackHeaders.h"
28#endif
29
30#include <math.h>
31#include <cstring>
32
33#if HAVE_LIBCOLPACK
34using namespace ColPack;
35#endif
36
37using namespace std;
38
39/****************************************************************************/
40/*******       sparse Jacobains, separate drivers             ***************/
41/****************************************************************************/
42
43/*--------------------------------------------------------------------------*/
44/*                                                sparsity pattern Jacobian */
45/*--------------------------------------------------------------------------*/
46/*                                                                         */
47
48int jac_pat(
49    short          tag,       /* tape identification                       */
50    int            depen,     /* number of dependent variables             */
51    int            indep,     /* number of independent variables           */
52    const double  *basepoint, /* independant variable values               */
53    unsigned int **crs,
54    /* returned compressed row block-index storage                         */
55    int          *options
56    /* control options
57                    options[0] : way of sparsity pattern computation
58                               0 - propagation of index domains (default)
59                               1 - propagation of bit pattern
60                    options[1] : test the computational graph control flow
61                               0 - safe mode (default)
62                               1 - tight mode
63                    options[2] : way of bit pattern propagation
64                               0 - automatic detection (default)
65                               1 - forward mode
66                               2 - reverse mode                            */
67) {
68    int             rc= -1;
69    int             i, ctrl_options[2];
70
71    if (crs == NULL) {
72        fprintf(DIAG_OUT,"ADOL-C user error in jac_pat(...) : "
73                "parameter crs may not be NULL !\n");
74        exit(-1);
75    } else
76        for (i=0; i<depen; i++)
77            crs[i] = NULL;
78
79    if (( options[0] < 0 ) || (options[0] > 1 ))
80        options[0] = 0; /* default */
81    if (( options[1] < 0 ) || (options[1] > 1 ))
82        options[1] = 0; /* default */
83    if (( options[2] < -1 ) || (options[2] > 2 ))
84        options[2] = 0; /* default */
85
86    if (options[0] == 0) {
87      if (options[1] == 1)
88        rc = indopro_forward_tight(tag, depen, indep, basepoint, crs);
89      else
90        {
91          rc = indopro_forward_safe(tag, depen, indep, basepoint, crs);
92        }
93    }
94    else 
95      {
96        ctrl_options[0] = options[1];
97        ctrl_options[1] = options[2];
98        rc = bit_vector_propagation( tag, depen, indep,  basepoint, crs, ctrl_options);
99      }
100
101    return(rc);
102}
103
104/*--------------------------------------------------------------------------*/
105/*                                                 seed matrix for Jacobian */
106/*--------------------------------------------------------------------------*/
107
108void generate_seed_jac
109(int m, int n, unsigned int **JP, double*** Seed, int *p, int option
110    /* control options
111                    option : way of compression
112                               0 - column compression (default)
113                               1 - row compression                */
114) 
115#if HAVE_LIBCOLPACK
116{
117  int dummy, i, j;
118
119  BipartiteGraphPartialColoringInterface *g = new BipartiteGraphPartialColoringInterface(SRC_MEM_ADOLC, JP, m, n);
120
121  if (option == 1) 
122    g->GenerateSeedJacobian_unmanaged(Seed, p, &dummy, 
123                                "SMALLEST_LAST","ROW_PARTIAL_DISTANCE_TWO"); 
124  else 
125    g->GenerateSeedJacobian_unmanaged(Seed, &dummy, p, 
126                                "SMALLEST_LAST","COLUMN_PARTIAL_DISTANCE_TWO"); 
127  delete g;
128
129}
130#else
131{
132    fprintf(DIAG_OUT, "ADOL-C error: function %s can only be used if linked with ColPack\n", __FUNCTION__);
133    exit(-1);
134}
135#endif
136
137/****************************************************************************/
138/*******        sparse Hessians, separate drivers             ***************/
139/****************************************************************************/
140
141/*---------------------------------------------------------------------------*/
142/*                                                  sparsity pattern Hessian */
143/*                                                                           */
144
145int hess_pat(
146    short          tag,        /* tape identification                        */
147    int            indep,      /* number of independent variables            */
148    const double  *basepoint,  /* independant variable values                */
149    unsigned int **crs,
150    /* returned compressed row block-index storage                         */
151    int          option
152    /* control option
153       option : test the computational graph control flow
154                               0 - safe mode (default)
155                               1 - tight mode                              */
156
157) {
158    int         rc= -1;
159    int         i;
160
161    if (crs == NULL) {
162        fprintf(DIAG_OUT,"ADOL-C user error in hess_pat(...) : "
163                "parameter crs may not be NULL !\n");
164        exit(-1);
165    } else
166        for (i=0; i<indep; i++)
167            crs[i] = NULL;
168
169    if (( option < 0 ) || (option > 2 ))
170      option = 0;   /* default */
171
172    if (option == 1)
173      rc = nonl_ind_forward_tight(tag, 1, indep, basepoint, crs);
174    else
175      rc = nonl_ind_forward_safe(tag, 1, indep, basepoint, crs);
176
177    return(rc);
178}
179
180/*--------------------------------------------------------------------------*/
181/*                                                  seed matrix for Hessian */
182/*--------------------------------------------------------------------------*/
183
184void generate_seed_hess
185(int n, unsigned int **HP, double*** Seed, int *p, int option
186    /* control options
187                    option : way of compression
188                               0 - indirect recovery (default)
189                               1 - direct recovery                */
190)
191#if HAVE_LIBCOLPACK
192{
193  int seed_rows, i, j;
194
195  GraphColoringInterface *g = new GraphColoringInterface(SRC_MEM_ADOLC, HP, n);
196
197  if (option == 0)
198    g->GenerateSeedHessian_unmanaged(Seed, &seed_rows, p, 
199                           "SMALLEST_LAST","ACYCLIC_FOR_INDIRECT_RECOVERY"); 
200  else
201    g->GenerateSeedHessian_unmanaged(Seed, &seed_rows, p, 
202                           "SMALLEST_LAST","STAR"); 
203  delete g;
204}
205#else
206{
207    fprintf(DIAG_OUT, "ADOL-C error: function %s can only be used if linked with ColPack\n", __FUNCTION__);
208    exit(-1);
209}
210#endif
211
212static void deepcopy_HP(unsigned int ***HPnew, unsigned int **HP, int indep)
213{
214    int i,j,s;
215    *HPnew = (unsigned int **)malloc(indep*sizeof(unsigned int *));
216    for (i=0; i<indep; i++) {
217       s=HP[i][0];
218       (*HPnew)[i] = (unsigned int *)malloc((s+1)*(sizeof(unsigned int)));
219       for (j=0; j<=s; j++)
220           (*HPnew)[i][j] = HP[i][j];
221    }
222}
223
224/****************************************************************************/
225/*******       sparse Jacobians, complete driver              ***************/
226/****************************************************************************/
227
228int sparse_jac(
229    short          tag,        /* tape identification                     */
230    int            depen,      /* number of dependent variables           */
231    int            indep,      /* number of independent variables         */
232    int            repeat,     /* indicated repeated call with same seed  */
233    const double  *basepoint,  /* independant variable values             */
234    int           *nnz,        /* number of nonzeros                      */
235    unsigned int **rind,       /* row index                               */
236    unsigned int **cind,       /* column index                            */
237    double       **values,     /* non-zero values                         */
238    int           *options
239    /* control options
240                    options[0] : way of sparsity pattern computation
241                               0 - propagation of index domains (default)
242                               1 - propagation of bit pattern
243                    options[1] : test the computational graph control flow
244                               0 - safe mode (default)
245                               1 - tight mode
246                    options[2] : way of bit pattern propagation
247                               0 - automatic detection (default)
248                               1 - forward mode
249                               2 - reverse mode                           
250                    options[3] : way of compression
251                               0 - column compression (default)
252                               1 - row compression                         */
253)
254#if HAVE_LIBCOLPACK
255{
256    int i;
257    unsigned int j;
258    SparseJacInfos sJinfos;
259    int dummy;
260    int ret_val;
261    BipartiteGraphPartialColoringInterface *g;
262    TapeInfos *tapeInfos;
263    JacobianRecovery1D *jr1d;
264    JacobianRecovery1D jr1d_loc;
265
266    ADOLC_OPENMP_THREAD_NUMBER;
267    ADOLC_OPENMP_GET_THREAD_NUMBER;
268
269    if (repeat == 0) {
270      if (( options[0] < 0 ) || (options[0] > 1 ))
271        options[0] = 0; /* default */
272      if (( options[1] < 0 ) || (options[1] > 1 ))
273        options[1] = 0; /* default */
274      if (( options[2] < -1 ) || (options[2] > 2 ))
275        options[2] = 0; /* default */
276      if (( options[3] < 0 ) || (options[3] > 1 ))
277        options[3] = 0; /* default */
278     
279      sJinfos.JP    = (unsigned int **) malloc(depen*sizeof(unsigned int *));
280      ret_val = jac_pat(tag, depen, indep, basepoint, sJinfos.JP, options);
281     
282
283      if (ret_val < 0) {
284        printf(" ADOL-C error in sparse_jac() \n");
285        return ret_val;
286      }
287     
288      sJinfos.depen = depen;
289      sJinfos.nnz_in = depen;
290      sJinfos.nnz_in = 0;
291      for (i=0;i<depen;i++) {
292            for (j=1;j<=sJinfos.JP[i][0];j++)
293              sJinfos.nnz_in++;
294      }
295     
296      *nnz = sJinfos.nnz_in;
297
298      if (options[2] == -1)
299        {
300          (*rind) = new unsigned int[*nnz];
301          (*cind) = new unsigned int[*nnz];
302          unsigned int index = 0;
303          for (i=0;i<depen;i++) 
304            for (j=1;j<=sJinfos.JP[i][0];j++)
305              {
306                (*rind)[index] = i;
307                (*cind)[index++] = sJinfos.JP[i][j];
308              }
309        }
310                       
311      /* sJinfos.Seed is memory managed by ColPack and will be deleted
312       * along with g. We only keep it in sJinfos for the repeat != 0 case */
313
314      g = new BipartiteGraphPartialColoringInterface(SRC_MEM_ADOLC, sJinfos.JP, depen, indep);
315      jr1d = new JacobianRecovery1D;
316       
317      if (options[3] == 1)
318        {
319          g->GenerateSeedJacobian(&(sJinfos.Seed), &(sJinfos.seed_rows), 
320                                  &(sJinfos.seed_clms), "SMALLEST_LAST","ROW_PARTIAL_DISTANCE_TWO"); 
321          sJinfos.seed_clms = indep;
322        } 
323      else
324        {
325          g->GenerateSeedJacobian(&(sJinfos.Seed), &(sJinfos.seed_rows), 
326                                &(sJinfos.seed_clms), "SMALLEST_LAST","COLUMN_PARTIAL_DISTANCE_TWO"); 
327          sJinfos.seed_rows = depen;
328        }
329     
330      sJinfos.B = myalloc2(sJinfos.seed_rows,sJinfos.seed_clms);
331      sJinfos.y = myalloc1(depen);
332     
333      sJinfos.g = (void *) g;
334      sJinfos.jr1d = (void *) jr1d;
335      setTapeInfoJacSparse(tag, sJinfos);
336      tapeInfos=getTapeInfos(tag);
337      memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
338    }
339    else
340      {
341        tapeInfos=getTapeInfos(tag);
342        memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
343        sJinfos.depen    = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.depen;
344        sJinfos.nnz_in    = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.nnz_in;
345        sJinfos.JP        = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.JP;
346        sJinfos.B         = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.B;
347        sJinfos.y         = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.y;
348        sJinfos.Seed      = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.Seed;
349        sJinfos.seed_rows = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.seed_rows;
350        sJinfos.seed_clms = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.seed_clms;
351        g = (BipartiteGraphPartialColoringInterface *)ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.g;
352        jr1d = (JacobianRecovery1D *)ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.jr1d;
353      }
354       
355    if (sJinfos.nnz_in != *nnz) {
356        printf(" ADOL-C error in sparse_jac():"
357               " Number of nonzeros not consistent,"
358               " repeat call with repeat = 0 \n");
359        return -3;
360    }
361
362    if (options[2] == -1)
363      return ret_val;
364
365    /* compute jacobian times matrix product */
366
367    if (options[3] == 1)
368      {
369        ret_val = zos_forward(tag,depen,indep,1,basepoint,sJinfos.y);
370        if (ret_val < 0) 
371          return ret_val;
372        MINDEC(ret_val,fov_reverse(tag,depen,indep,sJinfos.seed_rows,sJinfos.Seed,sJinfos.B));
373      }
374    else
375      ret_val = fov_forward(tag, depen, indep, sJinfos.seed_clms, basepoint, sJinfos.Seed, sJinfos.y, sJinfos.B);
376   
377
378    /* recover compressed Jacobian => ColPack library */
379 
380      if (options[3] == 1)
381       jr1d->RecoverD2Row_CoordinateFormat_unmanaged(g, sJinfos.B, sJinfos.JP, rind, cind, values);
382     else
383       jr1d->RecoverD2Cln_CoordinateFormat_unmanaged(g, sJinfos.B, sJinfos.JP, rind, cind, values);
384
385    return ret_val;
386
387}
388#else
389{
390    fprintf(DIAG_OUT, "ADOL-C error: function %s can only be used if linked with ColPack\n", __FUNCTION__);
391    exit(-1);
392}
393#endif
394
395/****************************************************************************/
396/*******        sparse Hessians, complete driver              ***************/
397/****************************************************************************/
398
399int sparse_hess(
400    short          tag,        /* tape identification                     */
401    int            indep,      /* number of independent variables         */
402    int            repeat,     /* indicated repeated call with same seed  */
403    const double  *basepoint,  /* independant variable values             */
404    int           *nnz,        /* number of nonzeros                      */
405    unsigned int **rind,       /* row index                               */
406    unsigned int **cind,       /* column index                            */
407    double       **values,     /* non-zero values                         */
408    int           *options
409    /* control options
410                    options[0] :test the computational graph control flow
411                               0 - safe mode (default)
412                               1 - tight mode
413                    options[1] : way of recovery
414                               0 - indirect recovery
415                               1 - direct recovery                         */
416)
417#if HAVE_LIBCOLPACK
418{
419    int i, l;
420    unsigned int j;
421    SparseHessInfos sHinfos;
422    double **Seed;
423    int dummy;
424    double y;
425    int ret_val=-1;
426    GraphColoringInterface *g;
427    TapeInfos *tapeInfos;
428    double *v, *w, **X, yt, lag=1;
429    HessianRecovery *hr;
430
431    ADOLC_OPENMP_THREAD_NUMBER;
432    ADOLC_OPENMP_GET_THREAD_NUMBER;
433
434    /* Generate sparsity pattern, determine nnz, allocate memory */
435    if (repeat <= 0) {
436        if (( options[0] < 0 ) || (options[0] > 1 ))
437          options[0] = 0; /* default */
438        if (( options[1] < 0 ) || (options[1] > 1 ))
439          options[1] = 0; /* default */
440
441        if (repeat == 0)
442          {
443            sHinfos.HP    = (unsigned int **) malloc(indep*sizeof(unsigned int *));
444
445            /* generate sparsity pattern */
446            ret_val = hess_pat(tag, indep, basepoint, sHinfos.HP, options[0]);
447
448            if (ret_val < 0) {
449              printf(" ADOL-C error in sparse_hess() \n");
450              return ret_val;
451            }
452          }
453        else
454          {
455            tapeInfos=getTapeInfos(tag);
456            memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
457            if (indep != ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.indep) {
458                fprintf(DIAG_OUT,"ADOL-C Error: wrong number of independents stored in hessian pattern.\n");
459                exit(-1);
460            }
461            deepcopy_HP(&sHinfos.HP,ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.HP,indep);         
462          }
463
464        sHinfos.indep = indep;
465        sHinfos.nnz_in = 0;
466
467        for (i=0;i<indep;i++) {
468            for (j=1;j<=sHinfos.HP[i][0];j++)
469              if ((int) sHinfos.HP[i][j] >= i)
470                sHinfos.nnz_in++;
471        }
472
473        *nnz = sHinfos.nnz_in;
474
475        /* compute seed matrix => ColPack library */
476
477        Seed = NULL;
478
479        g = new GraphColoringInterface(SRC_MEM_ADOLC, sHinfos.HP, indep);
480        hr = new HessianRecovery;
481
482        if (options[1] == 0)
483          g->GenerateSeedHessian(&Seed, &dummy, &sHinfos.p, 
484                                 "SMALLEST_LAST","ACYCLIC_FOR_INDIRECT_RECOVERY"); 
485        else
486          g->GenerateSeedHessian(&Seed, &dummy, &sHinfos.p, 
487                                 "SMALLEST_LAST","STAR"); 
488
489        sHinfos.Hcomp = myalloc2(indep,sHinfos.p);
490        sHinfos.Xppp = myalloc3(indep,sHinfos.p,1);
491
492        for (i=0; i<indep; i++)
493          for (l=0;l<sHinfos.p;l++)
494            sHinfos.Xppp[i][l][0] = Seed[i][l];
495
496        /* Seed will be freed by ColPack when g is freed */
497        Seed = NULL;
498
499        sHinfos.Yppp = myalloc3(1,sHinfos.p,1);
500
501        sHinfos.Zppp = myalloc3(sHinfos.p,indep,2);
502
503        sHinfos.Upp = myalloc2(1,2);
504        sHinfos.Upp[0][0] = 1;
505        sHinfos.Upp[0][1] = 0;
506
507        sHinfos.g = (void *) g;
508        sHinfos.hr = (void *) hr;
509
510        setTapeInfoHessSparse(tag, sHinfos);
511
512        tapeInfos=getTapeInfos(tag);
513        memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
514
515    }
516    else
517      {
518        tapeInfos=getTapeInfos(tag);
519        memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
520        sHinfos.nnz_in = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.nnz_in;
521        sHinfos.HP     = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.HP;
522        sHinfos.Hcomp  = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.Hcomp;
523        sHinfos.Xppp   = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.Xppp;
524        sHinfos.Yppp   = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.Yppp;
525        sHinfos.Zppp   = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.Zppp;
526        sHinfos.Upp    = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.Upp;
527        sHinfos.p      = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.p;
528        g = (GraphColoringInterface *)ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.g;
529        hr = (HessianRecovery *)ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.hr;
530      }
531
532    if (sHinfos.Upp == NULL) {
533        printf(" ADOL-C error in sparse_hess():"
534               " First call with repeat = 0 \n");
535        return -3;
536    }
537
538    if (sHinfos.nnz_in != *nnz) {
539        printf(" ADOL-C error in sparse_hess():"
540               " Number of nonzeros not consistent,"
541               " new call with repeat = 0 \n");
542        return -3;
543    }
544
545    if (repeat == -1)
546      return ret_val;
547
548//     this is the most efficient variant. However, there is somewhere a bug in hos_ov_reverse
549//     ret_val = hov_wk_forward(tag,1,indep,1,2,sHinfos.p,basepoint,sHinfos.Xppp,&y,sHinfos.Yppp);
550//     MINDEC(ret_val,hos_ov_reverse(tag,1,indep,1,sHinfos.p,sHinfos.Upp,sHinfos.Zppp));
551
552//     for (i = 0; i < sHinfos.p; ++i)
553//       for (l = 0; l < indep; ++l)
554//      sHinfos.Hcomp[l][i] = sHinfos.Zppp[i][l][1];
555
556//
557//     therefore, we use hess_vec isntead of hess_mat
558
559    v    = (double*) malloc(indep*sizeof(double));
560    w    = (double*) malloc(indep*sizeof(double));
561    X = myalloc2(indep,2);
562//         sHinfos.Xppp = myalloc3(indep,sHinfos.p,1);
563
564    for (i = 0; i < sHinfos.p; ++i)
565      {
566        for (l = 0; l < indep; ++l)
567          {
568            v[l] = sHinfos.Xppp[l][i][0];
569          }
570        ret_val = fos_forward(tag, 1, indep, 2, basepoint, v, &y, &yt);
571        MINDEC(ret_val, hos_reverse(tag, 1, indep, 1, &lag, X));
572        for (l = 0; l < indep; ++l)
573          {
574            sHinfos.Hcomp[l][i] = X[l][1];
575          }
576      }
577
578    myfree1(v);
579    myfree1(w);
580    myfree2(X);   
581
582
583    /* recover compressed Hessian => ColPack library */
584
585//      if (options[1] == 0)
586//        HessianRecovery::IndirectRecover_CoordinateFormat(g, sHinfos.Hcomp, sHinfos.HP, rind, cind, values);
587//      else
588//        HessianRecovery::DirectRecover_CoordinateFormat(g, sHinfos.Hcomp, sHinfos.HP, rind, cind, values);
589 
590     if (options[1] == 0)
591       hr->IndirectRecover_CoordinateFormat_unmanaged(g, sHinfos.Hcomp, sHinfos.HP, rind, cind, values);
592     else
593       hr->DirectRecover_CoordinateFormat_unmanaged(g, sHinfos.Hcomp, sHinfos.HP, rind, cind, values);
594 
595    return ret_val;
596
597}
598#else
599{
600    fprintf(DIAG_OUT, "ADOL-C error: function %s can only be used if linked with ColPack\n", __FUNCTION__);
601    exit(-1);
602}
603#endif
604
605
606/****************************************************************************/
607/*******      sparse Hessians, set and get sparsity pattern   ***************/
608/****************************************************************************/
609
610void set_HP(
611    short          tag,        /* tape identification                     */
612    int            indep,      /* number of independent variables         */
613    unsigned int ** HP)
614#ifdef SPARSE
615{
616    SparseHessInfos sHinfos;
617    TapeInfos *tapeInfos;
618
619    tapeInfos=getTapeInfos(tag);
620    memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
621    sHinfos.nnz_in = 0;
622    deepcopy_HP(&sHinfos.HP,HP,indep);
623    sHinfos.Hcomp  = NULL;
624    sHinfos.Xppp   = NULL;
625    sHinfos.Yppp   = NULL;
626    sHinfos.Zppp   = NULL;
627    sHinfos.Upp    = NULL;
628    sHinfos.p      = 0;
629    sHinfos.g      = NULL;
630    sHinfos.hr     = NULL;
631    sHinfos.indep  = indep;
632    setTapeInfoHessSparse(tag, sHinfos);
633}
634#else
635{
636    fprintf(DIAG_OUT, "ADOL-C error: function %s can only be used if sparse configuration option was used\n", __FUNCTION__);
637    exit(-1);
638}
639#endif
640
641void get_HP(
642    short          tag,        /* tape identification                     */
643    int            indep,      /* number of independent variables         */
644    unsigned int *** HP)
645#ifdef SPARSE
646{
647    SparseHessInfos sHinfos;
648    TapeInfos *tapeInfos;
649
650    tapeInfos=getTapeInfos(tag);
651    memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
652    deepcopy_HP(HP,ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.HP,indep);
653}
654#else
655{
656    fprintf(DIAG_OUT, "ADOL-C error: function %s can only be used if sparse configuration option was used\n", __FUNCTION__);
657    exit(-1);
658}
659#endif
660
661/*****************************************************************************/
662/*                                                    JACOBIAN BLOCK PATTERN */
663
664/* ------------------------------------------------------------------------- */
665int bit_vector_propagation(
666    short          tag,        /* tape identification                */
667    int            depen,      /* number of dependent variables      */
668    int            indep,      /* number of independent variables    */
669    const double  *basepoint, /* independant variable values         */
670    unsigned int **crs,
671    /* compressed block row storage                                  */
672    int *options       /* control options                            */
673    /* options[0] : way of bit pattern propagation
674                    0 - automatic detection (default)
675                    1 - forward mode
676                    2 - reverse mode   
677       options[1] : test the computational graph control flow
678                    0 - safe variant (default)
679                    1 - tight variant  */   
680) {
681
682    int                rc= 3;
683    char               forward_mode, tight_mode;
684    int                i, ii, j, jj, k, k_old, bits_per_long,
685    i_blocks_per_strip, d_blocks_per_strip;
686    int                this_strip_i_bl_idx, next_strip_i_bl_idx,
687    this_strip_d_bl_idx, next_strip_d_bl_idx;
688    int                stripmined_calls, strip_idx;
689    int                p_stripmine, q_stripmine, p_ind_bl_bp, q_dep_bl_bp,
690    i_bl_idx, d_bl_idx;
691    unsigned long int  value1, v;
692    unsigned long int  **seed=NULL, *s, **jac_bit_pat=NULL, *jac;
693    unsigned char      *indep_blocks_flags=NULL, *i_b_flags;
694    double             *valuepoint=NULL;
695
696   if ( options[1] == 0 ) {
697        if ( depen >= indep/2 )
698            options[1] = 1; /* forward */
699        else
700            options[1] = 2; /* reverse */
701    }
702
703    if ( options[1] == 1 )
704        forward_mode = 1;
705    else
706        forward_mode = 0;
707
708    if ( options[0] == 1 )
709        tight_mode = 1;
710    else
711        tight_mode = 0;
712
713    if ( ! forward_mode )
714        valuepoint = myalloc1(depen);
715
716    /* bit pattern parameters */
717
718    /* number of bits in an unsigned long int variable */
719    bits_per_long = 8 * sizeof(unsigned long int);
720    /* olvo 20000214 nl: inserted explicit cast to unsigned long int */
721    value1 =  (unsigned long int) 1 << (bits_per_long - 1); /* 10000....0 */
722
723    /* =================================================== forward propagation */
724    if ( forward_mode ) {
725     
726        if (( tight_mode ) && ( basepoint == NULL )) {
727            fprintf(DIAG_OUT, "ADOL-C error in jac_pat(...) :  supply basepoint x for tight mode.\n");
728            exit(-1);
729        }
730
731        /* indep partial derivatives for the whole Jacobian */
732
733        /* number of unsigned longs to store the whole seed / Jacobian matrice */
734        p_ind_bl_bp = indep / bits_per_long
735                      + ( (indep % bits_per_long) != 0 );
736
737        /* number of unsigned longs to store the seed / Jacobian strips */
738        if ( p_ind_bl_bp <= PQ_STRIPMINE_MAX ) {
739            p_stripmine = p_ind_bl_bp;
740            stripmined_calls = 1;
741        } else {
742            p_stripmine = PQ_STRIPMINE_MAX;
743            stripmined_calls = p_ind_bl_bp / PQ_STRIPMINE_MAX
744                               + ( (p_ind_bl_bp % PQ_STRIPMINE_MAX) != 0 );
745        }
746
747        /* number of independent blocks per seed / Jacobian strip */
748        i_blocks_per_strip = p_stripmine * bits_per_long;
749
750        /* allocate memory --------------------------------------------------- */
751
752        if ( ! (indep_blocks_flags = (unsigned char*)
753                                     calloc(i_blocks_per_strip, sizeof(char)) ) ) {
754            fprintf(DIAG_OUT, "ADOL-C error, "__FILE__
755                    ":%i : \njac_pat(...) unable to allocate %i bytes !\n",
756                    __LINE__, (int)(i_blocks_per_strip*sizeof(char)));
757            exit(-1);
758        }
759
760        seed        = myalloc2_ulong(indep, p_stripmine);
761        jac_bit_pat = myalloc2_ulong(depen, p_stripmine);
762
763        /* strip-mining : repeated forward calls ----------------------------- */
764
765        for (strip_idx = 0; strip_idx < stripmined_calls; strip_idx++) {
766            /* build a partition of the seed matrix (indep x indep_blocks) --- */
767            /* (indep x i_blocks_per_strip) as a bit pattern                   */
768            s = seed[0];
769            for (i=0; i<indep; i++)
770                for (ii=0; ii<p_stripmine; ii++) /* 2 loops if short -> int !!! */
771                    *s++ = 0; /* set old seed matrix to 0 */
772
773            this_strip_i_bl_idx = strip_idx * i_blocks_per_strip;
774            next_strip_i_bl_idx = (strip_idx+1) * i_blocks_per_strip;
775            if ( next_strip_i_bl_idx > indep )
776                next_strip_i_bl_idx = indep;
777            v = value1; /* 10000....0 */
778
779            for (i=0; i<indep; i++)
780                if ( (this_strip_i_bl_idx <= i)
781                        && (i < next_strip_i_bl_idx) ) {
782                    ii = (i - this_strip_i_bl_idx) /  bits_per_long;
783                    seed[i][ii] = v >> ((i - this_strip_i_bl_idx) %  bits_per_long);
784                }
785
786            /* bit pattern propagation by forward ---------------------------- */
787
788            if ( tight_mode )
789              {
790                rc = int_forward_tight( tag, depen, indep, p_stripmine,
791                                        basepoint, seed, valuepoint, jac_bit_pat);
792              }
793            else
794              {
795                rc = int_forward_safe ( tag, depen, indep, p_stripmine,
796                                        seed, jac_bit_pat);
797              }
798
799            /* extract  pattern from bit patterns --------------------- */
800
801            for (j = 0; j < depen; j++) {
802                    ii = -1;
803                    v = 0;
804
805                    jac = jac_bit_pat[j];
806                    i_b_flags = indep_blocks_flags;
807                    for (i_bl_idx = 0; i_bl_idx < i_blocks_per_strip; i_bl_idx++) {
808                        if ( !v ) {
809                            v =  value1; /* 10000....0 */
810                            ii++;
811                        }
812                        if ( v & jac[ii] )
813                            *i_b_flags = 1;
814                        i_b_flags++;
815
816                        v = v >> 1;
817                    }
818
819                if ( strip_idx == 0 )
820                    k_old = 0;
821                else
822                    k_old = crs[j][0];
823                k = 0;
824                i_b_flags = indep_blocks_flags;
825                for (i = 0; i < i_blocks_per_strip; i++)
826                    k += *i_b_flags++;
827
828                if ((k > 0 ) || ( strip_idx == 0 )) {
829                    if ( ! (crs[j] = (unsigned int*)realloc(crs[j],
830                                            (k_old+k+1)*sizeof(unsigned int))) ) {
831                        fprintf(DIAG_OUT, "ADOL-C error, "__FILE__
832                                 ":%i : \njac_pat(...) unable to allocate %i bytes !\n",
833                                __LINE__, (int)((k_old+k+1)*sizeof(unsigned int)));
834                        exit(-1);
835                    }
836                    if ( strip_idx == 0 )
837                        crs[j][0]  = 0;
838                    if ( k > 0 ) {
839                        k = crs[j][0] + 1;
840                        i_b_flags = indep_blocks_flags;
841                        for (i = 0; i < i_blocks_per_strip; i++) {
842                            if ( *i_b_flags ) {
843                                crs[j][k++] = this_strip_i_bl_idx + i;
844                                *i_b_flags = 0;
845                            }
846                            i_b_flags++;
847                        }
848                        /* current/total number of non-zero blocks of indep. vars. */
849                        crs[j][0] = k - 1;
850                    }
851                }
852            }
853        } /* strip_idx */
854
855    } /* forward */
856
857
858    /* =================================================== reverse propagation */
859    else {
860
861        /* depen weight vectors for the whole Jacobian */
862
863        /* number of unsigned longs to store the whole seed / Jacobian matrice */
864        q_dep_bl_bp = depen / bits_per_long
865                      + ( (depen % bits_per_long) != 0 );
866
867        /* number of unsigned longs to store the seed / Jacobian strips */
868        if ( q_dep_bl_bp <= PQ_STRIPMINE_MAX ) {
869            q_stripmine = q_dep_bl_bp;
870            stripmined_calls = 1;
871        } else {
872            q_stripmine = PQ_STRIPMINE_MAX;
873            stripmined_calls = q_dep_bl_bp / PQ_STRIPMINE_MAX
874                               + ( (q_dep_bl_bp % PQ_STRIPMINE_MAX) != 0 );
875        }
876
877        /* number of dependent blocks per seed / Jacobian strip */
878        d_blocks_per_strip = q_stripmine * bits_per_long;
879
880        /* allocate memory --------------------------------------------------- */
881        if ( ! (indep_blocks_flags = (unsigned char*)calloc(indep,
882                                     sizeof(unsigned char)) ) ) {
883            fprintf(DIAG_OUT, "ADOL-C error, "__FILE__
884                    ":%i : \njac_pat(...) unable to allocate %i bytes !\n",
885                    __LINE__, (int)(indep*sizeof(unsigned char)));
886            exit(-1);
887        }
888
889        seed        = myalloc2_ulong(q_stripmine, depen);
890        jac_bit_pat = myalloc2_ulong(q_stripmine, indep);
891
892
893        /* olvo 20000214: call to forward required in tight mode only,
894           in safe mode no basepoint available! */
895        if ( tight_mode ) {
896            if ( basepoint == NULL ) {
897                fprintf(DIAG_OUT, "ADOL-C error in jac_pat(..) :  ");
898                fprintf(DIAG_OUT, "no basepoint x for tight mode supplied.\n");
899                exit(-1);
900            }
901
902            rc = zos_forward(tag, depen, indep, 1, basepoint, valuepoint);
903        }
904
905        /* strip-mining : repeated reverse calls ----------------------------- */
906
907        for (strip_idx = 0; strip_idx < stripmined_calls; strip_idx++) {
908            /* build a partition of the seed matrix (depen_blocks x depen)     */
909            /* (d_blocks_per_strip x depen) as a bit pattern                   */
910            s = seed[0];
911            for (jj=0; jj<q_stripmine; jj++) /* 2 loops if short -> int !!! */
912                for (j=0; j<depen; j++)
913                    *s++ = 0; /* set old seed matrix to 0 */
914
915            this_strip_d_bl_idx = strip_idx * d_blocks_per_strip;
916            next_strip_d_bl_idx = (strip_idx+1) * d_blocks_per_strip;
917            if ( next_strip_d_bl_idx > depen )
918                next_strip_d_bl_idx = depen;
919            v = value1; /* 10000....0 */
920
921            for (j=0; j<depen; j++)
922                if ( (this_strip_d_bl_idx <= j)
923                        && (j < next_strip_d_bl_idx) ) {
924                    jj = (j - this_strip_d_bl_idx) /  bits_per_long;
925                    seed[jj][j] = v >> ((j - this_strip_d_bl_idx) % bits_per_long);
926                }
927
928            /* bit pattern propagation by reverse ---------------------------- */
929
930            if ( tight_mode )
931                rc = int_reverse_tight( tag, depen, indep, q_stripmine,
932                                        seed, jac_bit_pat);
933            else
934                rc = int_reverse_safe ( tag, depen, indep, q_stripmine,
935                                        seed, jac_bit_pat);
936
937
938            /* extract pattern from bit patterns --------------------- */
939
940            jj = -1;
941            v = 0;
942            for (d_bl_idx = this_strip_d_bl_idx;
943                    d_bl_idx < next_strip_d_bl_idx; d_bl_idx++) {
944                if ( !v ) {
945                    v =  value1; /* 10000....0 */
946                    jj++;
947                }
948                jac = jac_bit_pat[jj];
949                for (i=0; i<indep; i++) {
950                    if ( v & *jac++ ) {
951                        indep_blocks_flags[i] = 1;
952                    }
953                }
954
955                v = v >> 1;
956
957                k=0;
958                i_b_flags = indep_blocks_flags;
959                for (i=0; i<indep; i++)
960                    k += *i_b_flags++;
961
962                if ( ! (crs[d_bl_idx] = (unsigned int*)malloc((k+1)*sizeof(unsigned int))) ) {
963                    fprintf(DIAG_OUT, "ADOL-C error, "__FILE__
964                            ":%i : \njac_pat(...) unable to allocate %i bytes !\n",
965                            __LINE__, (int)((k+1)*sizeof(unsigned int)));
966                    exit(-1);
967                }
968                crs[d_bl_idx][0] = k; /* number of non-zero indep. blocks */
969                k=1;
970                i_b_flags = indep_blocks_flags;
971                for (i=0; i<indep; i++) {
972                    if ( *i_b_flags ) {
973                        crs[d_bl_idx][k++] = i;
974                        *i_b_flags = 0;
975                    }
976                    i_b_flags++;
977                }
978            }
979
980        } /* strip_idx */
981
982    } /* reverse */
983
984    if ( ! forward_mode ) {
985        free((char*)valuepoint);
986        valuepoint=NULL;
987    }
988    free((char*)*seed);
989    free((char*)seed);
990    seed=NULL;
991    free((char*)*jac_bit_pat);
992    free((char*)jac_bit_pat);
993    jac_bit_pat=NULL;
994    free((char*)indep_blocks_flags);
995    indep_blocks_flags=NULL;
996
997    return(rc);
998}
999
1000BEGIN_C_DECLS
1001/*****************************************************************************/
1002/*                                                FREE SPARSE JACOBIAN INFOS */
1003
1004/* ------------------------------------------------------------------------- */
1005void freeSparseJacInfos(double *y, double **B, unsigned int **JP, void *g, 
1006                        void *jr1d, int seed_rows, int seed_clms, int depen)
1007{
1008    int i;
1009    if(y)
1010      myfree1(y);
1011
1012    if(B)
1013      myfree2(B);
1014
1015    for (int i=0;i<depen;i++) {
1016      free(JP[i]);
1017    }
1018
1019    free(JP);
1020
1021#ifdef HAVE_LIBCOLPACK
1022     if (g) 
1023       delete (BipartiteGraphPartialColoringInterface *) g;
1024
1025    if (jr1d)
1026        delete (JacobianRecovery1D*)jr1d;
1027#endif
1028}
1029/*****************************************************************************/
1030/*                                                 FREE SPARSE HESSIAN INFOS */
1031
1032/* ------------------------------------------------------------------------- */
1033void freeSparseHessInfos(double **Hcomp, double ***Xppp, double ***Yppp, double ***Zppp, 
1034                         double **Upp, unsigned int **HP,
1035                         void *g, void *hr, int p, int indep)
1036{
1037    int i;
1038
1039    if(Hcomp)
1040      myfree2(Hcomp);
1041
1042   if(Xppp)
1043      myfree3(Xppp);
1044
1045   if(Yppp)
1046      myfree3(Yppp);
1047
1048   if(Zppp)
1049      myfree3(Zppp);
1050   if(Upp)
1051      myfree2(Upp);
1052
1053   if(HP)
1054     {
1055       for (int i=0;i<indep;i++) {
1056         free(HP[i]);
1057       }
1058       free(HP);
1059     }
1060
1061#ifdef HAVE_LIBCOLPACK
1062     if (g) 
1063       delete (GraphColoringInterface *) g;
1064    if (hr)
1065        delete (HessianRecovery*) hr;
1066#endif
1067}
1068
1069END_C_DECLS
1070
1071/****************************************************************************/
1072/*                                                               THAT'S ALL */
1073
Note: See TracBrowser for help on using the repository browser.