source: trunk/ADOL-C/examples/additional_examples/sparse/sparse_jacobian.cpp @ 171

Last change on this file since 171 was 171, checked in by kulshres, 9 years ago

Squashed merge branch 'master' of 'gitclone' into svn

  • 'master' of 'gitclone': (84 commits) adjust example makefiles and include paths get rid of the symlink in the src subdirectory

details of the commits:
commit c9e4bc332d2363f737fc2e8a8fcfc2e43ddb9d15
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon Oct 4 15:43:47 2010 +0200

adjust example makefiles and include paths

include paths in example sources were wrong for some time now
simplify makefile rules too, there is really no need for checking SPARSE
adjust include paths in makefiles.

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

commit e6e1963e41e097fd5b4a79cd1611c12f6868dc94
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon Oct 4 15:41:25 2010 +0200

get rid of the symlink in the src subdirectory

windows doesn't like symlinks and make infinite depth directories
we now create a symlink for build in the directory parallel to src
adjust all makefiles.am accordingly for build

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

  • Property svn:keywords set to Author Date Id Revision
File size: 8.5 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     sparse_jacobian.cpp
4 Revision: $Id: sparse_jacobian.cpp 171 2010-10-04 13:57:19Z kulshres $
5 Contents: example for computation of sparse jacobians
6
7 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
8               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel
9 
10 This file is part of ADOL-C. This software is provided as open source.
11 Any use, reproduction, or distribution of the software constitutes
12 recipient's acceptance of the terms of the accompanying license file.
13 
14---------------------------------------------------------------------------*/
15
16#include <math.h>
17#include <cstdlib>
18#include <cstdio>
19
20#include <adolc/adolc.h>
21#include <adolc/adolc_sparse.h>
22
23#define tag 1
24
25void   ceval_ad(adouble *x, adouble *c);
26void   ceval(double *x, double *c);
27
28void printmat(char* name, int n, int m, double** M);
29
30int main() {
31    int n=6, m=3;
32    double x[6], c[3];
33    adouble xad[6], cad[3];
34
35    int i, j;
36
37/****************************************************************************/
38/*******                function evaluation                   ***************/
39/****************************************************************************/
40
41    for(i=0;i<n;i++)
42        x[i] = log(1.0+i);
43
44    /* Tracing of function c(x) */
45
46    trace_on(tag);
47      for(i=0;i<n;i++)
48        xad[i] <<= x[i];
49
50      ceval_ad(xad,cad);
51
52      for(i=0;i<m;i++)
53        cad[i] >>= c[i];
54    trace_off();
55
56    printf("\n c =  ");
57    for(j=0;j<m;j++)
58        printf(" %e ",c[j]);
59    printf("\n");
60
61/****************************************************************************/
62/********           For comparisons: Full Jacobian                   ********/
63/****************************************************************************/
64
65    double **J;
66    J = myalloc2(m,n);
67
68    jacobian(tag,m,n,x,J);
69
70    printmat(" J",m,n,J);
71    printf("\n");
72
73
74/****************************************************************************/
75/*******       sparse Jacobians, complete driver              ***************/
76/****************************************************************************/
77
78    /* coordinate format for Jacobian */
79    unsigned int *rind  = NULL;        /* row indices    */
80    unsigned int *cind  = NULL;        /* column indices */
81    double       *values = NULL;       /* values         */
82    int nnz;
83    int options[4];
84
85    options[0] = 0;          /* sparsity pattern by index domains (default) */ 
86    options[1] = 0;          /*                         safe mode (default) */ 
87    options[2] = 0;          /*              not required if options[0] = 0 */ 
88    options[3] = 0;          /*                column compression (default) */ 
89
90    sparse_jac(tag, m, n, 0, x, &nnz, &rind, &cind, &values, options);
91
92    printf("In sparse format:\n");
93    for (i=0;i<nnz;i++)
94        printf("%2d %2d %10.6f\n\n",rind[i],cind[i],values[i]);
95
96    delete[] rind;
97    delete[] cind;
98    delete[] values;
99/*--------------------------------------------------------------------------*/
100/*  same approach but using row compression                                 */
101/*--------------------------------------------------------------------------*/
102
103    options[3] = 1;                   /*   row compression => reverse mode, */ 
104                                      /* sometimes better than forward mode */ 
105                                      /* due to sparsity structure          */
106
107    sparse_jac(tag, m, n, 0, x, &nnz, &rind, &cind, &values, options);
108
109    printf("In sparse format (using row compression): \n");
110    for (i=0;i<nnz;i++)
111        printf("%2d %2d %10.6f\n\n",rind[i],cind[i],values[i]);
112
113    delete[] rind;
114    delete[] cind;
115    delete[] values;
116/*--------------------------------------------------------------------------*/
117/*  change value of x, but not the sparsity pattern                         */
118/*--------------------------------------------------------------------------*/
119
120    for(i=0;i<n;i++)
121        x[i] = 2.0*i;
122
123/*  For comparisons: Full Jacobian                                          */
124
125    jacobian(tag,m,n,x,J);
126
127    printmat(" J",m,n,J);
128    printf("\n");
129
130/*  repeated call of sparse_jac with same sparsity pattern => repeat = 1 */
131
132    sparse_jac(tag, m, n, 1, x, &nnz, &rind, &cind, &values, options);
133
134    printf("In sparse format:\n");
135    for (i=0;i<nnz;i++)
136        printf("%2d %2d %10.6f\n\n",rind[i],cind[i],values[i]);
137
138    delete[] rind;
139    delete[] cind;
140    delete[] values;
141/*--------------------------------------------------------------------------*/
142/*  same approach but using row compression                                 */
143/*--------------------------------------------------------------------------*/
144
145    options[3] = 1;                   /*   row compression => reverse mode, */ 
146                                      /* sometimes better than forward mode */ 
147                                      /* due to sparsity structure          */
148
149    sparse_jac(tag, m, n, 0, x, &nnz, &rind, &cind, &values, options);
150
151    printf("In sparse format (using row compression): \n");
152    for (i=0;i<nnz;i++)
153        printf("%2d %2d %10.6f\n\n",rind[i],cind[i],values[i]);
154
155    delete[] rind;
156    delete[] cind;
157    delete[] values;
158/****************************************************************************/
159/*******       sparse Jacobians, separate drivers             ***************/
160/****************************************************************************/
161
162/*--------------------------------------------------------------------------*/
163/*                                                sparsity pattern Jacobian */
164/*--------------------------------------------------------------------------*/
165
166    unsigned int  **JP=NULL;                /* compressed block row storage */
167    int ctrl[3];
168
169    JP = (unsigned int **) malloc(m*sizeof(unsigned int*));
170    ctrl[0] = 0;
171    ctrl[1] = 0;
172    ctrl[2] = 0;
173
174    jac_pat(tag, m, n, x, JP, ctrl);
175
176    printf("\n");
177    printf("Sparsity pattern of Jacobian: \n");
178    for (i=0;i<m;i++) {
179        printf(" %d: ",i);
180        for (j=1;j<= (int) JP[i][0];j++)
181            printf(" %d ",JP[i][j]);
182        printf("\n");
183    }
184    printf("\n");
185
186
187/*--------------------------------------------------------------------------*/
188/*                                                              seed matrix */
189/*--------------------------------------------------------------------------*/
190
191    double **Seed;
192    int p;
193    int option = 0;
194
195    /* option = 0 column compression (default),
196       option = 1 rom compression                */
197 
198    generate_seed_jac(m, n, JP, &Seed, &p, option);
199
200    printf(" p_J = %d \n",p);
201    printmat(" Seed matrix",n,p,Seed);
202    printf("\n");
203
204/*--------------------------------------------------------------------------*/
205/*                                                      compressed Jacobian */
206/*--------------------------------------------------------------------------*/
207
208    double **Jcomp;
209    Jcomp = myalloc2(m,p);
210
211    fov_forward(tag,m,n,p,x,Seed,c,Jcomp);
212    printmat("compressed J:",m,p,Jcomp);
213    printf("\n");
214
215
216/*--------------------------------------------------------------------------*/
217/*  change value of x, but not the sparsity pattern                         */
218/*--------------------------------------------------------------------------*/
219
220    for(i=0;i<n;i++)
221        x[i] = 2.0*i;
222
223/*  For comparisons: Full Jacobian                                          */
224
225    jacobian(tag,m,n,x,J);
226
227    printmat(" J",m,n,J);
228    printf("\n");
229
230
231    fov_forward(tag,m,n,p,x,Seed,c,Jcomp);
232    printmat("compressed J:",m,p,Jcomp);
233    printf("\n");
234
235    for (i=0;i<m;i++)
236        free(JP[i]);
237    free(JP);
238    myfree2(J);
239
240    for (i = 0; i < n; i++)
241        delete[] Seed[i];
242    delete[] Seed;
243
244    myfree2(Jcomp);
245}
246
247
248/***************************************************************************/
249
250void ceval(double *x, double *c) {
251    c[0] = 2*x[0]+x[1]-2.0;
252    c[1] = x[2]*x[2]+x[3]*x[3]-2.0;
253    c[2] = 3*x[4]*x[5] - 3.0;
254}
255
256/***************************************************************************/
257
258void ceval_ad(adouble *x, adouble *c) {
259    c[0] = 2*x[0]+x[1]-2.0;
260    c[0] += cos(x[3])*sin(x[4]);
261    c[1] = x[2]*x[2]+x[3]*x[3]-2.0;
262    c[2] = 3*x[4]*x[5] - 3.0+sin(x[4]*x[5]);
263}
264
265/***************************************************************************/
266
267void printmat(char* name, int m, int n, double** M) {
268    int i,j;
269
270    printf("%s \n",name);
271    for(i=0; i<m ;i++) {
272        printf("\n %d: ",i);
273        for(j=0;j<n ;j++)
274            printf(" %10.4f ", M[i][j]);
275    }
276    printf("\n");
277}
Note: See TracBrowser for help on using the repository browser.