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

Last change on this file since 324 was 324, checked in by awalther, 8 years ago

delete timings in sparsedrivers.cpp

  • Property svn:keywords set to Author Date Id Revision
File size: 38.4 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     sparse/sparsedrivers.cpp
4 Revision: $Id: sparsedrivers.cpp 324 2012-06-01 14:19:18Z 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 <ColPack/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) = (unsigned int*)calloc(*nnz,sizeof(unsigned int));
301          (*cind) = (unsigned int*)calloc(*nnz,sizeof(unsigned int));
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_OMP"); 
321          sJinfos.seed_clms = indep;
322          ret_val = sJinfos.seed_rows;
323        } 
324      else
325        {
326          g->GenerateSeedJacobian(&(sJinfos.Seed), &(sJinfos.seed_rows), 
327                                &(sJinfos.seed_clms), "SMALLEST_LAST","COLUMN_PARTIAL_DISTANCE_TWO_OMP"); 
328          sJinfos.seed_rows = depen;
329          ret_val = sJinfos.seed_clms;
330        }
331     
332      sJinfos.B = myalloc2(sJinfos.seed_rows,sJinfos.seed_clms);
333      sJinfos.y = myalloc1(depen);
334     
335      sJinfos.g = (void *) g;
336      sJinfos.jr1d = (void *) jr1d;
337      setTapeInfoJacSparse(tag, sJinfos);
338      tapeInfos=getTapeInfos(tag);
339      memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
340    }
341    else
342      {
343        tapeInfos=getTapeInfos(tag);
344        memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
345        sJinfos.depen    = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.depen;
346        sJinfos.nnz_in    = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.nnz_in;
347        sJinfos.JP        = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.JP;
348        sJinfos.B         = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.B;
349        sJinfos.y         = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.y;
350        sJinfos.Seed      = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.Seed;
351        sJinfos.seed_rows = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.seed_rows;
352        sJinfos.seed_clms = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.seed_clms;
353        g = (BipartiteGraphPartialColoringInterface *)ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.g;
354        jr1d = (JacobianRecovery1D *)ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sJinfos.jr1d;
355      }
356       
357    if (sJinfos.nnz_in != *nnz) {
358        printf(" ADOL-C error in sparse_jac():"
359               " Number of nonzeros not consistent,"
360               " repeat call with repeat = 0 \n");
361        return -3;
362    }
363
364    if (options[2] == -1)
365      return ret_val;
366
367    /* compute jacobian times matrix product */
368
369    if (options[3] == 1)
370      {
371        ret_val = zos_forward(tag,depen,indep,1,basepoint,sJinfos.y);
372        if (ret_val < 0) 
373          return ret_val;
374        MINDEC(ret_val,fov_reverse(tag,depen,indep,sJinfos.seed_rows,sJinfos.Seed,sJinfos.B));
375      }
376    else
377      ret_val = fov_forward(tag, depen, indep, sJinfos.seed_clms, basepoint, sJinfos.Seed, sJinfos.y, sJinfos.B);
378   
379
380    /* recover compressed Jacobian => ColPack library */
381 
382    if (*values != NULL && *rind != NULL && *cind != NULL) {
383      // everything is preallocated, we assume correctly
384      // call usermem versions
385      if (options[3] == 1)
386       jr1d->RecoverD2Row_CoordinateFormat_usermem_OMP(g, sJinfos.B, sJinfos.JP, rind, cind, values);
387     else
388       jr1d->RecoverD2Cln_CoordinateFormat_usermem_OMP(g, sJinfos.B, sJinfos.JP, rind, cind, values);
389    } else {
390      // at least one of rind cind values is not allocated, deallocate others
391      // and call unmanaged versions
392      if (*values != NULL)
393          free(*values);
394      if (*rind != NULL)
395          free(*rind);
396      if (*cind != NULL)
397          free(*cind);
398      if (options[3] == 1)
399       jr1d->RecoverD2Row_CoordinateFormat_unmanaged_OMP(g, sJinfos.B, sJinfos.JP, rind, cind, values);
400     else
401       jr1d->RecoverD2Cln_CoordinateFormat_unmanaged_OMP(g, sJinfos.B, sJinfos.JP, rind, cind, values);
402    }
403
404    return ret_val;
405
406}
407#else
408{
409    fprintf(DIAG_OUT, "ADOL-C error: function %s can only be used if linked with ColPack\n", __FUNCTION__);
410    exit(-1);
411}
412#endif
413
414/****************************************************************************/
415/*******        sparse Hessians, complete driver              ***************/
416/****************************************************************************/
417
418int sparse_hess(
419    short          tag,        /* tape identification                     */
420    int            indep,      /* number of independent variables         */
421    int            repeat,     /* indicated repeated call with same seed  */
422    const double  *basepoint,  /* independant variable values             */
423    int           *nnz,        /* number of nonzeros                      */
424    unsigned int **rind,       /* row index                               */
425    unsigned int **cind,       /* column index                            */
426    double       **values,     /* non-zero values                         */
427    int           *options
428    /* control options
429                    options[0] :test the computational graph control flow
430                               0 - safe mode (default)
431                               1 - tight mode
432                    options[1] : way of recovery
433                               0 - indirect recovery
434                               1 - direct recovery                         */
435)
436#if HAVE_LIBCOLPACK
437{
438    int i, l;
439    unsigned int j;
440    SparseHessInfos sHinfos;
441    double **Seed;
442    int dummy;
443    double y;
444    int ret_val=-1;
445    GraphColoringInterface *g;
446    TapeInfos *tapeInfos;
447    double *v, *w, **X, yt, lag=1;
448    HessianRecovery *hr;
449
450    ADOLC_OPENMP_THREAD_NUMBER;
451    ADOLC_OPENMP_GET_THREAD_NUMBER;
452
453    /* Generate sparsity pattern, determine nnz, allocate memory */
454    if (repeat <= 0) {
455        if (( options[0] < 0 ) || (options[0] > 1 ))
456          options[0] = 0; /* default */
457        if (( options[1] < 0 ) || (options[1] > 1 ))
458          options[1] = 0; /* default */
459
460        if (repeat == 0)
461          {
462            sHinfos.HP    = (unsigned int **) malloc(indep*sizeof(unsigned int *));
463
464            /* generate sparsity pattern */
465            ret_val = hess_pat(tag, indep, basepoint, sHinfos.HP, options[0]);
466
467            if (ret_val < 0) {
468              printf(" ADOL-C error in sparse_hess() \n");
469              return ret_val;
470            }
471          }
472        else
473          {
474            tapeInfos=getTapeInfos(tag);
475            memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
476            if (indep != ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.indep) {
477                fprintf(DIAG_OUT,"ADOL-C Error: wrong number of independents stored in hessian pattern.\n");
478                exit(-1);
479            }
480            deepcopy_HP(&sHinfos.HP,ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.HP,indep);         
481          }
482
483        sHinfos.indep = indep;
484        sHinfos.nnz_in = 0;
485
486        for (i=0;i<indep;i++) {
487            for (j=1;j<=sHinfos.HP[i][0];j++)
488              if ((int) sHinfos.HP[i][j] >= i)
489                sHinfos.nnz_in++;
490        }
491
492        *nnz = sHinfos.nnz_in;
493
494        /* compute seed matrix => ColPack library */
495
496        Seed = NULL;
497
498        g = new GraphColoringInterface(SRC_MEM_ADOLC, sHinfos.HP, indep);
499        hr = new HessianRecovery;
500
501        if (options[1] == 0)
502          g->GenerateSeedHessian(&Seed, &dummy, &sHinfos.p, 
503                                 "SMALLEST_LAST","ACYCLIC_FOR_INDIRECT_RECOVERY"); 
504        else
505          g->GenerateSeedHessian(&Seed, &dummy, &sHinfos.p, 
506                                 "SMALLEST_LAST","STAR"); 
507
508        sHinfos.Hcomp = myalloc2(indep,sHinfos.p);
509        sHinfos.Xppp = myalloc3(indep,sHinfos.p,1);
510
511        for (i=0; i<indep; i++)
512          for (l=0;l<sHinfos.p;l++)
513            sHinfos.Xppp[i][l][0] = Seed[i][l];
514
515        /* Seed will be freed by ColPack when g is freed */
516        Seed = NULL;
517
518        sHinfos.Yppp = myalloc3(1,sHinfos.p,1);
519
520        sHinfos.Zppp = myalloc3(sHinfos.p,indep,2);
521
522        sHinfos.Upp = myalloc2(1,2);
523        sHinfos.Upp[0][0] = 1;
524        sHinfos.Upp[0][1] = 0;
525
526        sHinfos.g = (void *) g;
527        sHinfos.hr = (void *) hr;
528
529        setTapeInfoHessSparse(tag, sHinfos);
530
531        tapeInfos=getTapeInfos(tag);
532        memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
533
534    }
535    else
536      {
537        tapeInfos=getTapeInfos(tag);
538        memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
539        sHinfos.nnz_in = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.nnz_in;
540        sHinfos.HP     = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.HP;
541        sHinfos.Hcomp  = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.Hcomp;
542        sHinfos.Xppp   = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.Xppp;
543        sHinfos.Yppp   = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.Yppp;
544        sHinfos.Zppp   = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.Zppp;
545        sHinfos.Upp    = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.Upp;
546        sHinfos.p      = ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.p;
547        g = (GraphColoringInterface *)ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.g;
548        hr = (HessianRecovery *)ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.hr;
549      }
550
551    if (sHinfos.Upp == NULL) {
552        printf(" ADOL-C error in sparse_hess():"
553               " First call with repeat = 0 \n");
554        return -3;
555    }
556
557    if (sHinfos.nnz_in != *nnz) {
558        printf(" ADOL-C error in sparse_hess():"
559               " Number of nonzeros not consistent,"
560               " new call with repeat = 0 \n");
561        return -3;
562    }
563
564    if (repeat == -1)
565      return ret_val;
566
567//     this is the most efficient variant. However, there is somewhere a bug in hos_ov_reverse
568//     ret_val = hov_wk_forward(tag,1,indep,1,2,sHinfos.p,basepoint,sHinfos.Xppp,&y,sHinfos.Yppp);
569//     MINDEC(ret_val,hos_ov_reverse(tag,1,indep,1,sHinfos.p,sHinfos.Upp,sHinfos.Zppp));
570
571//     for (i = 0; i < sHinfos.p; ++i)
572//       for (l = 0; l < indep; ++l)
573//      sHinfos.Hcomp[l][i] = sHinfos.Zppp[i][l][1];
574
575//
576//     therefore, we use hess_vec isntead of hess_mat
577
578    v    = (double*) malloc(indep*sizeof(double));
579    w    = (double*) malloc(indep*sizeof(double));
580    X = myalloc2(indep,2);
581//         sHinfos.Xppp = myalloc3(indep,sHinfos.p,1);
582
583    for (i = 0; i < sHinfos.p; ++i)
584      {
585        for (l = 0; l < indep; ++l)
586          {
587            v[l] = sHinfos.Xppp[l][i][0];
588          }
589        ret_val = fos_forward(tag, 1, indep, 2, basepoint, v, &y, &yt);
590        MINDEC(ret_val, hos_reverse(tag, 1, indep, 1, &lag, X));
591        for (l = 0; l < indep; ++l)
592          {
593            sHinfos.Hcomp[l][i] = X[l][1];
594          }
595      }
596
597    myfree1(v);
598    myfree1(w);
599    myfree2(X);   
600
601
602    /* recover compressed Hessian => ColPack library */
603
604//      if (options[1] == 0)
605//        HessianRecovery::IndirectRecover_CoordinateFormat(g, sHinfos.Hcomp, sHinfos.HP, rind, cind, values);
606//      else
607//        HessianRecovery::DirectRecover_CoordinateFormat(g, sHinfos.Hcomp, sHinfos.HP, rind, cind, values);
608 
609    if (*values != NULL && *rind != NULL && *cind != NULL) {
610     // everything is preallocated, we assume correctly
611     // call usermem versions
612     if (options[1] == 0)
613       hr->IndirectRecover_CoordinateFormat_usermem(g, sHinfos.Hcomp, sHinfos.HP, rind, cind, values);
614     else
615       hr->DirectRecover_CoordinateFormat_usermem_OMP(g, sHinfos.Hcomp, sHinfos.HP, rind, cind, values);
616    } else {
617      // at least one of rind cind values is not allocated, deallocate others
618      // and call unmanaged versions
619      if (*values != NULL)
620          free(*values);
621      if (*rind != NULL)
622          free(*rind);
623      if (*cind != NULL)
624          free(*cind);
625     if (options[1] == 0)
626       hr->IndirectRecover_CoordinateFormat_unmanaged(g, sHinfos.Hcomp, sHinfos.HP, rind, cind, values);
627     else
628       hr->DirectRecover_CoordinateFormat_unmanaged_OMP(g, sHinfos.Hcomp, sHinfos.HP, rind, cind, values);
629    }
630    return ret_val;
631
632}
633#else
634{
635    fprintf(DIAG_OUT, "ADOL-C error: function %s can only be used if linked with ColPack\n", __FUNCTION__);
636    exit(-1);
637}
638#endif
639
640
641/****************************************************************************/
642/*******      sparse Hessians, set and get sparsity pattern   ***************/
643/****************************************************************************/
644
645void set_HP(
646    short          tag,        /* tape identification                     */
647    int            indep,      /* number of independent variables         */
648    unsigned int ** HP)
649#ifdef SPARSE
650{
651    SparseHessInfos sHinfos;
652    TapeInfos *tapeInfos;
653
654    ADOLC_OPENMP_THREAD_NUMBER;
655    ADOLC_OPENMP_GET_THREAD_NUMBER;
656
657    tapeInfos=getTapeInfos(tag);
658    memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
659    sHinfos.nnz_in = 0;
660    deepcopy_HP(&sHinfos.HP,HP,indep);
661    sHinfos.Hcomp  = NULL;
662    sHinfos.Xppp   = NULL;
663    sHinfos.Yppp   = NULL;
664    sHinfos.Zppp   = NULL;
665    sHinfos.Upp    = NULL;
666    sHinfos.p      = 0;
667    sHinfos.g      = NULL;
668    sHinfos.hr     = NULL;
669    sHinfos.indep  = indep;
670    setTapeInfoHessSparse(tag, sHinfos);
671}
672#else
673{
674    fprintf(DIAG_OUT, "ADOL-C error: function %s can only be used if sparse configuration option was used\n", __FUNCTION__);
675    exit(-1);
676}
677#endif
678
679void get_HP(
680    short          tag,        /* tape identification                     */
681    int            indep,      /* number of independent variables         */
682    unsigned int *** HP)
683#ifdef SPARSE
684{
685    SparseHessInfos sHinfos;
686    TapeInfos *tapeInfos;
687
688    ADOLC_OPENMP_THREAD_NUMBER;
689    ADOLC_OPENMP_GET_THREAD_NUMBER;
690
691    tapeInfos=getTapeInfos(tag);
692    memcpy(&ADOLC_CURRENT_TAPE_INFOS, tapeInfos, sizeof(TapeInfos));
693    deepcopy_HP(HP,ADOLC_CURRENT_TAPE_INFOS.pTapeInfos.sHinfos.HP,indep);
694}
695#else
696{
697    fprintf(DIAG_OUT, "ADOL-C error: function %s can only be used if sparse configuration option was used\n", __FUNCTION__);
698    exit(-1);
699}
700#endif
701
702/*****************************************************************************/
703/*                                                    JACOBIAN BLOCK PATTERN */
704
705/* ------------------------------------------------------------------------- */
706int bit_vector_propagation(
707    short          tag,        /* tape identification                */
708    int            depen,      /* number of dependent variables      */
709    int            indep,      /* number of independent variables    */
710    const double  *basepoint, /* independant variable values         */
711    unsigned int **crs,
712    /* compressed block row storage                                  */
713    int *options       /* control options                            */
714    /* options[0] : way of bit pattern propagation
715                    0 - automatic detection (default)
716                    1 - forward mode
717                    2 - reverse mode   
718       options[1] : test the computational graph control flow
719                    0 - safe variant (default)
720                    1 - tight variant  */   
721) {
722
723    int                rc= 3;
724    char               forward_mode, tight_mode;
725    int                i, ii, j, jj, k, k_old, bits_per_long,
726    i_blocks_per_strip, d_blocks_per_strip;
727    int                this_strip_i_bl_idx, next_strip_i_bl_idx,
728    this_strip_d_bl_idx, next_strip_d_bl_idx;
729    int                stripmined_calls, strip_idx;
730    int                p_stripmine, q_stripmine, p_ind_bl_bp, q_dep_bl_bp,
731    i_bl_idx, d_bl_idx;
732    unsigned long int  value1, v;
733    unsigned long int  **seed=NULL, *s, **jac_bit_pat=NULL, *jac;
734    unsigned char      *indep_blocks_flags=NULL, *i_b_flags;
735    double             *valuepoint=NULL;
736
737   if ( options[1] == 0 ) {
738        if ( depen >= indep/2 )
739            options[1] = 1; /* forward */
740        else
741            options[1] = 2; /* reverse */
742    }
743
744    if ( options[1] == 1 )
745        forward_mode = 1;
746    else
747        forward_mode = 0;
748
749    if ( options[0] == 1 )
750        tight_mode = 1;
751    else
752        tight_mode = 0;
753
754    if ( ! forward_mode )
755        valuepoint = myalloc1(depen);
756
757    /* bit pattern parameters */
758
759    /* number of bits in an unsigned long int variable */
760    bits_per_long = 8 * sizeof(unsigned long int);
761    /* olvo 20000214 nl: inserted explicit cast to unsigned long int */
762    value1 =  (unsigned long int) 1 << (bits_per_long - 1); /* 10000....0 */
763
764    /* =================================================== forward propagation */
765    if ( forward_mode ) {
766     
767        if (( tight_mode ) && ( basepoint == NULL )) {
768            fprintf(DIAG_OUT, "ADOL-C error in jac_pat(...) :  supply basepoint x for tight mode.\n");
769            exit(-1);
770        }
771
772        /* indep partial derivatives for the whole Jacobian */
773
774        /* number of unsigned longs to store the whole seed / Jacobian matrice */
775        p_ind_bl_bp = indep / bits_per_long
776                      + ( (indep % bits_per_long) != 0 );
777
778        /* number of unsigned longs to store the seed / Jacobian strips */
779        if ( p_ind_bl_bp <= PQ_STRIPMINE_MAX ) {
780            p_stripmine = p_ind_bl_bp;
781            stripmined_calls = 1;
782        } else {
783            p_stripmine = PQ_STRIPMINE_MAX;
784            stripmined_calls = p_ind_bl_bp / PQ_STRIPMINE_MAX
785                               + ( (p_ind_bl_bp % PQ_STRIPMINE_MAX) != 0 );
786        }
787
788        /* number of independent blocks per seed / Jacobian strip */
789        i_blocks_per_strip = p_stripmine * bits_per_long;
790
791        /* allocate memory --------------------------------------------------- */
792
793        if ( ! (indep_blocks_flags = (unsigned char*)
794                                     calloc(i_blocks_per_strip, sizeof(char)) ) ) {
795            fprintf(DIAG_OUT, "ADOL-C error, "__FILE__
796                    ":%i : \njac_pat(...) unable to allocate %i bytes !\n",
797                    __LINE__, (int)(i_blocks_per_strip*sizeof(char)));
798            exit(-1);
799        }
800
801        seed        = myalloc2_ulong(indep, p_stripmine);
802        jac_bit_pat = myalloc2_ulong(depen, p_stripmine);
803
804        /* strip-mining : repeated forward calls ----------------------------- */
805
806        for (strip_idx = 0; strip_idx < stripmined_calls; strip_idx++) {
807            /* build a partition of the seed matrix (indep x indep_blocks) --- */
808            /* (indep x i_blocks_per_strip) as a bit pattern                   */
809            s = seed[0];
810            for (i=0; i<indep; i++)
811                for (ii=0; ii<p_stripmine; ii++) /* 2 loops if short -> int !!! */
812                    *s++ = 0; /* set old seed matrix to 0 */
813
814            this_strip_i_bl_idx = strip_idx * i_blocks_per_strip;
815            next_strip_i_bl_idx = (strip_idx+1) * i_blocks_per_strip;
816            if ( next_strip_i_bl_idx > indep )
817                next_strip_i_bl_idx = indep;
818            v = value1; /* 10000....0 */
819
820            for (i=0; i<indep; i++)
821                if ( (this_strip_i_bl_idx <= i)
822                        && (i < next_strip_i_bl_idx) ) {
823                    ii = (i - this_strip_i_bl_idx) /  bits_per_long;
824                    seed[i][ii] = v >> ((i - this_strip_i_bl_idx) %  bits_per_long);
825                }
826
827            /* bit pattern propagation by forward ---------------------------- */
828
829            if ( tight_mode )
830              {
831                rc = int_forward_tight( tag, depen, indep, p_stripmine,
832                                        basepoint, seed, valuepoint, jac_bit_pat);
833              }
834            else
835              {
836                rc = int_forward_safe ( tag, depen, indep, p_stripmine,
837                                        seed, jac_bit_pat);
838              }
839
840            /* extract  pattern from bit patterns --------------------- */
841
842            for (j = 0; j < depen; j++) {
843                    ii = -1;
844                    v = 0;
845
846                    jac = jac_bit_pat[j];
847                    i_b_flags = indep_blocks_flags;
848                    for (i_bl_idx = 0; i_bl_idx < i_blocks_per_strip; i_bl_idx++) {
849                        if ( !v ) {
850                            v =  value1; /* 10000....0 */
851                            ii++;
852                        }
853                        if ( v & jac[ii] )
854                            *i_b_flags = 1;
855                        i_b_flags++;
856
857                        v = v >> 1;
858                    }
859
860                if ( strip_idx == 0 )
861                    k_old = 0;
862                else
863                    k_old = crs[j][0];
864                k = 0;
865                i_b_flags = indep_blocks_flags;
866                for (i = 0; i < i_blocks_per_strip; i++)
867                    k += *i_b_flags++;
868
869                if ((k > 0 ) || ( strip_idx == 0 )) {
870                    if ( ! (crs[j] = (unsigned int*)realloc(crs[j],
871                                            (k_old+k+1)*sizeof(unsigned int))) ) {
872                        fprintf(DIAG_OUT, "ADOL-C error, "__FILE__
873                                 ":%i : \njac_pat(...) unable to allocate %i bytes !\n",
874                                __LINE__, (int)((k_old+k+1)*sizeof(unsigned int)));
875                        exit(-1);
876                    }
877                    if ( strip_idx == 0 )
878                        crs[j][0]  = 0;
879                    if ( k > 0 ) {
880                        k = crs[j][0] + 1;
881                        i_b_flags = indep_blocks_flags;
882                        for (i = 0; i < i_blocks_per_strip; i++) {
883                            if ( *i_b_flags ) {
884                                crs[j][k++] = this_strip_i_bl_idx + i;
885                                *i_b_flags = 0;
886                            }
887                            i_b_flags++;
888                        }
889                        /* current/total number of non-zero blocks of indep. vars. */
890                        crs[j][0] = k - 1;
891                    }
892                }
893            }
894        } /* strip_idx */
895
896    } /* forward */
897
898
899    /* =================================================== reverse propagation */
900    else {
901
902        /* depen weight vectors for the whole Jacobian */
903
904        /* number of unsigned longs to store the whole seed / Jacobian matrice */
905        q_dep_bl_bp = depen / bits_per_long
906                      + ( (depen % bits_per_long) != 0 );
907
908        /* number of unsigned longs to store the seed / Jacobian strips */
909        if ( q_dep_bl_bp <= PQ_STRIPMINE_MAX ) {
910            q_stripmine = q_dep_bl_bp;
911            stripmined_calls = 1;
912        } else {
913            q_stripmine = PQ_STRIPMINE_MAX;
914            stripmined_calls = q_dep_bl_bp / PQ_STRIPMINE_MAX
915                               + ( (q_dep_bl_bp % PQ_STRIPMINE_MAX) != 0 );
916        }
917
918        /* number of dependent blocks per seed / Jacobian strip */
919        d_blocks_per_strip = q_stripmine * bits_per_long;
920
921        /* allocate memory --------------------------------------------------- */
922        if ( ! (indep_blocks_flags = (unsigned char*)calloc(indep,
923                                     sizeof(unsigned char)) ) ) {
924            fprintf(DIAG_OUT, "ADOL-C error, "__FILE__
925                    ":%i : \njac_pat(...) unable to allocate %i bytes !\n",
926                    __LINE__, (int)(indep*sizeof(unsigned char)));
927            exit(-1);
928        }
929
930        seed        = myalloc2_ulong(q_stripmine, depen);
931        jac_bit_pat = myalloc2_ulong(q_stripmine, indep);
932
933
934        /* olvo 20000214: call to forward required in tight mode only,
935           in safe mode no basepoint available! */
936        if ( tight_mode ) {
937            if ( basepoint == NULL ) {
938                fprintf(DIAG_OUT, "ADOL-C error in jac_pat(..) :  ");
939                fprintf(DIAG_OUT, "no basepoint x for tight mode supplied.\n");
940                exit(-1);
941            }
942
943            rc = zos_forward(tag, depen, indep, 1, basepoint, valuepoint);
944        }
945
946        /* strip-mining : repeated reverse calls ----------------------------- */
947
948        for (strip_idx = 0; strip_idx < stripmined_calls; strip_idx++) {
949            /* build a partition of the seed matrix (depen_blocks x depen)     */
950            /* (d_blocks_per_strip x depen) as a bit pattern                   */
951            s = seed[0];
952            for (jj=0; jj<q_stripmine; jj++) /* 2 loops if short -> int !!! */
953                for (j=0; j<depen; j++)
954                    *s++ = 0; /* set old seed matrix to 0 */
955
956            this_strip_d_bl_idx = strip_idx * d_blocks_per_strip;
957            next_strip_d_bl_idx = (strip_idx+1) * d_blocks_per_strip;
958            if ( next_strip_d_bl_idx > depen )
959                next_strip_d_bl_idx = depen;
960            v = value1; /* 10000....0 */
961
962            for (j=0; j<depen; j++)
963                if ( (this_strip_d_bl_idx <= j)
964                        && (j < next_strip_d_bl_idx) ) {
965                    jj = (j - this_strip_d_bl_idx) /  bits_per_long;
966                    seed[jj][j] = v >> ((j - this_strip_d_bl_idx) % bits_per_long);
967                }
968
969            /* bit pattern propagation by reverse ---------------------------- */
970
971            if ( tight_mode )
972                rc = int_reverse_tight( tag, depen, indep, q_stripmine,
973                                        seed, jac_bit_pat);
974            else
975                rc = int_reverse_safe ( tag, depen, indep, q_stripmine,
976                                        seed, jac_bit_pat);
977
978
979            /* extract pattern from bit patterns --------------------- */
980
981            jj = -1;
982            v = 0;
983            for (d_bl_idx = this_strip_d_bl_idx;
984                    d_bl_idx < next_strip_d_bl_idx; d_bl_idx++) {
985                if ( !v ) {
986                    v =  value1; /* 10000....0 */
987                    jj++;
988                }
989                jac = jac_bit_pat[jj];
990                for (i=0; i<indep; i++) {
991                    if ( v & *jac++ ) {
992                        indep_blocks_flags[i] = 1;
993                    }
994                }
995
996                v = v >> 1;
997
998                k=0;
999                i_b_flags = indep_blocks_flags;
1000                for (i=0; i<indep; i++)
1001                    k += *i_b_flags++;
1002
1003                if ( ! (crs[d_bl_idx] = (unsigned int*)malloc((k+1)*sizeof(unsigned int))) ) {
1004                    fprintf(DIAG_OUT, "ADOL-C error, "__FILE__
1005                            ":%i : \njac_pat(...) unable to allocate %i bytes !\n",
1006                            __LINE__, (int)((k+1)*sizeof(unsigned int)));
1007                    exit(-1);
1008                }
1009                crs[d_bl_idx][0] = k; /* number of non-zero indep. blocks */
1010                k=1;
1011                i_b_flags = indep_blocks_flags;
1012                for (i=0; i<indep; i++) {
1013                    if ( *i_b_flags ) {
1014                        crs[d_bl_idx][k++] = i;
1015                        *i_b_flags = 0;
1016                    }
1017                    i_b_flags++;
1018                }
1019            }
1020
1021        } /* strip_idx */
1022
1023    } /* reverse */
1024
1025    if ( ! forward_mode ) {
1026        free((char*)valuepoint);
1027        valuepoint=NULL;
1028    }
1029    free((char*)*seed);
1030    free((char*)seed);
1031    seed=NULL;
1032    free((char*)*jac_bit_pat);
1033    free((char*)jac_bit_pat);
1034    jac_bit_pat=NULL;
1035    free((char*)indep_blocks_flags);
1036    indep_blocks_flags=NULL;
1037
1038    return(rc);
1039}
1040
1041BEGIN_C_DECLS
1042/*****************************************************************************/
1043/*                                                FREE SPARSE JACOBIAN INFOS */
1044
1045/* ------------------------------------------------------------------------- */
1046void freeSparseJacInfos(double *y, double **B, unsigned int **JP, void *g, 
1047                        void *jr1d, int seed_rows, int seed_clms, int depen)
1048{
1049    int i;
1050    if(y)
1051      myfree1(y);
1052
1053    if(B)
1054      myfree2(B);
1055
1056    for (int i=0;i<depen;i++) {
1057      free(JP[i]);
1058    }
1059
1060    free(JP);
1061
1062#ifdef HAVE_LIBCOLPACK
1063     if (g) 
1064       delete (BipartiteGraphPartialColoringInterface *) g;
1065
1066    if (jr1d)
1067        delete (JacobianRecovery1D*)jr1d;
1068#endif
1069}
1070/*****************************************************************************/
1071/*                                                 FREE SPARSE HESSIAN INFOS */
1072
1073/* ------------------------------------------------------------------------- */
1074void freeSparseHessInfos(double **Hcomp, double ***Xppp, double ***Yppp, double ***Zppp, 
1075                         double **Upp, unsigned int **HP,
1076                         void *g, void *hr, int p, int indep)
1077{
1078    int i;
1079
1080    if(Hcomp)
1081      myfree2(Hcomp);
1082
1083   if(Xppp)
1084      myfree3(Xppp);
1085
1086   if(Yppp)
1087      myfree3(Yppp);
1088
1089   if(Zppp)
1090      myfree3(Zppp);
1091
1092   if(Upp)
1093      myfree2(Upp);
1094
1095   if(HP)
1096     {
1097       for (int i=0;i<indep;i++) 
1098         free(HP[i]);
1099       free(HP);
1100     }
1101
1102#ifdef HAVE_LIBCOLPACK
1103     if (g) 
1104       delete (GraphColoringInterface *) g;
1105    if (hr)
1106        delete (HessianRecovery*) hr;
1107#endif
1108}
1109
1110END_C_DECLS
1111
1112/****************************************************************************/
1113/*                                                               THAT'S ALL */
1114
Note: See TracBrowser for help on using the repository browser.