[40] | 1 | /*---------------------------------------------------------------------------- |
---|
| 2 | ADOL-C -- Automatic Differentiation by Overloading in C++ |
---|
| 3 | File: LUdet.cpp |
---|
[42] | 4 | Revision: $Id: LUdet.cpp 88 2010-02-17 18:33:52Z utke $ |
---|
[40] | 5 | Contents: example for |
---|
| 6 | * Computation of the determinant of a matrix |
---|
| 7 | by LU-decomposition of the system matrix without pivoting |
---|
| 8 | * application of tapedoc to observe taping of |
---|
| 9 | the new op_codes for the elementary operations |
---|
| 10 | |
---|
| 11 | y += x1 * x2; |
---|
| 12 | y -= x1 * x2; |
---|
| 13 | |
---|
| 14 | Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz, |
---|
| 15 | Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel |
---|
| 16 | |
---|
| 17 | This file is part of ADOL-C. This software is provided as open source. |
---|
| 18 | Any use, reproduction, or distribution of the software constitutes |
---|
| 19 | recipient's acceptance of the terms of the accompanying license file. |
---|
| 20 | |
---|
| 21 | ---------------------------------------------------------------------------*/ |
---|
| 22 | |
---|
| 23 | /****************************************************************************/ |
---|
| 24 | /* INCLUDES */ |
---|
| 25 | #include <../examples/additional_examples/lufact/LU.h> |
---|
| 26 | |
---|
| 27 | |
---|
| 28 | /****************************************************************************/ |
---|
| 29 | /* MAIN PROGRAM */ |
---|
| 30 | int main() { /*------------------------------------------------------------------------*/ |
---|
| 31 | /* variables */ |
---|
| 32 | const int tag = 1; // tape tag |
---|
| 33 | const int size = 5; // system size |
---|
| 34 | const int indep = size*size; // # of indeps |
---|
| 35 | const int depen = 1; // # of deps |
---|
| 36 | |
---|
| 37 | double A[size][size], a1[size], a2[size], det; // passive variables |
---|
| 38 | adouble **AA, *AAp, Adet; // active variables |
---|
| 39 | double *args = myalloc1(indep); // arguments |
---|
| 40 | double *grad = myalloc1(indep); // the gradient |
---|
| 41 | double **hess = myalloc2(indep,indep); // the hessian |
---|
| 42 | |
---|
| 43 | int i,j; |
---|
| 44 | |
---|
| 45 | |
---|
| 46 | /*------------------------------------------------------------------------*/ |
---|
| 47 | /* Info */ |
---|
| 48 | fprintf(stdout,"DETERMINANT by LU-DECOMPOSITION (ADOL-C Example)\n\n"); |
---|
| 49 | |
---|
| 50 | |
---|
| 51 | /*------------------------------------------------------------------------*/ |
---|
| 52 | /* Allcoation und initialization of the system matrix */ |
---|
| 53 | AA = new adouble*[size]; |
---|
| 54 | AAp = new adouble[size*size]; |
---|
| 55 | for (i=0; i<size; i++) { |
---|
| 56 | AA[i] = AAp; |
---|
| 57 | AAp += size; |
---|
| 58 | } |
---|
| 59 | for(i=0; i<size; i++) { |
---|
| 60 | a1[i] = i*0.25; |
---|
| 61 | a2[i] = i*0.33; |
---|
| 62 | } |
---|
| 63 | for(i=0; i<size; i++) { |
---|
| 64 | for(j=0; j<size; j++) |
---|
| 65 | A[i][j] = a1[i]*a2[j]; |
---|
| 66 | A[i][i] += i+1; |
---|
| 67 | } |
---|
| 68 | |
---|
| 69 | |
---|
| 70 | /*------------------------------------------------------------------------*/ |
---|
| 71 | /* Taping the computation of the determinant */ |
---|
| 72 | trace_on(tag); |
---|
| 73 | /* marking indeps */ |
---|
| 74 | for(i=0; i<size; i++) |
---|
| 75 | for(j=0; j<size; j++) |
---|
| 76 | AA[i][j] <<= (args[i*size+j] = A[i][j]); // indep. vars |
---|
| 77 | /* LU-factorization and computation of determinant */ |
---|
| 78 | LUfact(size,AA); |
---|
| 79 | Adet = AA[0][0]; |
---|
| 80 | for (i=1; i<size; i++) |
---|
| 81 | Adet *= AA[i][i]; |
---|
| 82 | /* marking deps */ |
---|
| 83 | Adet >>= det; |
---|
| 84 | trace_off(); |
---|
[88] | 85 | fprintf(stdout," Determinant (original): %16.4E\n",det); |
---|
[40] | 86 | |
---|
| 87 | |
---|
| 88 | /*------------------------------------------------------------------------*/ |
---|
| 89 | /* Recomputation of determinant */ |
---|
| 90 | function(tag,depen,indep,args,&det); |
---|
[88] | 91 | fprintf(stdout," Determinant (from tape): %16.4E\n",det); |
---|
[40] | 92 | |
---|
| 93 | |
---|
| 94 | /*------------------------------------------------------------------------*/ |
---|
| 95 | /* Computation of gradient */ |
---|
| 96 | gradient(tag,indep,args,grad); |
---|
| 97 | fprintf(stdout," Gradient:\n"); |
---|
| 98 | for (i=0; i<size; i++) { |
---|
| 99 | for (j=0; j<size; j++) |
---|
[88] | 100 | fprintf(stdout," %14.6E",grad[i*size+j]); |
---|
[40] | 101 | fprintf(stdout,"\n"); |
---|
| 102 | } |
---|
| 103 | |
---|
| 104 | /*------------------------------------------------------------------------*/ |
---|
| 105 | /* Computation of hessian */ |
---|
| 106 | hessian(tag,indep,args,hess); |
---|
| 107 | fprintf(stdout," Part of Hessian:\n"); |
---|
| 108 | for (i=0; i<size; i++) { |
---|
| 109 | for (j=0; j<size; j++) |
---|
[88] | 110 | fprintf(stdout," %14.6E",hess[0][i*size+j]); |
---|
[40] | 111 | fprintf(stdout,"\n"); |
---|
| 112 | } |
---|
| 113 | |
---|
| 114 | /*------------------------------------------------------------------------*/ |
---|
| 115 | /* Tape-documentation */ |
---|
| 116 | tape_doc(tag,depen,indep,args,&det); |
---|
| 117 | |
---|
| 118 | /*------------------------------------------------------------------------*/ |
---|
| 119 | /* Tape statistics */ |
---|
| 120 | int tape_stats[STAT_SIZE]; |
---|
| 121 | tapestats(tag,tape_stats); |
---|
| 122 | |
---|
| 123 | fprintf(stdout,"\n independents %d\n",tape_stats[NUM_INDEPENDENTS]); |
---|
| 124 | fprintf(stdout," dependents %d\n",tape_stats[NUM_DEPENDENTS]); |
---|
| 125 | fprintf(stdout," operations %d\n",tape_stats[NUM_OPERATIONS]); |
---|
| 126 | fprintf(stdout," operations buffer size %d\n",tape_stats[OP_BUFFER_SIZE]); |
---|
| 127 | fprintf(stdout," locations buffer size %d\n",tape_stats[LOC_BUFFER_SIZE]); |
---|
| 128 | fprintf(stdout," constants buffer size %d\n",tape_stats[VAL_BUFFER_SIZE]); |
---|
| 129 | fprintf(stdout," maxlive %d\n",tape_stats[NUM_MAX_LIVES]); |
---|
| 130 | fprintf(stdout," valstack size %d\n\n",tape_stats[TAY_STACK_SIZE]); |
---|
| 131 | |
---|
| 132 | /*------------------------------------------------------------------------*/ |
---|
| 133 | /* That's it */ |
---|
| 134 | return 1; |
---|
| 135 | } |
---|
| 136 | |
---|
| 137 | |
---|
| 138 | |
---|
| 139 | |
---|
| 140 | |
---|
| 141 | |
---|
| 142 | |
---|
| 143 | |
---|