1 2 #ifndef lint 3 static char vcid[] = "$Id: fdmatrix.c,v 1.2 1996/11/07 15:09:08 bsmith Exp bsmith $"; 4 #endif 5 6 /* 7 This is where the abstract matrix operations are defined that are 8 used for finite difference computations of Jacobians using coloring. 9 */ 10 11 #include "petsc.h" 12 #include "src/mat/matimpl.h" /*I "mat.h" I*/ 13 #include "src/vec/vecimpl.h" 14 #include "pinclude/pviewer.h" 15 16 /*@C 17 MatFDColoringView - Views a finite difference coloring context. 18 19 Input Parameter: 20 . color - the coloring context 21 22 23 .seealso: MatFDColoringCreate() 24 @*/ 25 int MatFDColoringView(MatFDColoring color,Viewer viewer) 26 { 27 int i,j,format,ierr; 28 29 if (!viewer) viewer = VIEWER_STDOUT_WORLD; 30 31 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 32 33 if (format == VIEWER_FORMAT_ASCII_INFO) { 34 printf("MatFDColoring Object:\n"); 35 printf(" Error tolerance %g\n",color->error_rel); 36 printf(" umin %g\n",color->umin); 37 } else { 38 printf("MatFDColoring Object:\n"); 39 printf(" Error tolerance %g\n",color->error_rel); 40 printf(" umin %g\n",color->umin); 41 printf("Number of colors %d\n",color->ncolors); 42 for ( i=0; i<color->ncolors; i++ ) { 43 printf("Information for color %d\n",i); 44 printf("Number of columns %d\n",color->ncolumns[i]); 45 for ( j=0; j<color->ncolumns[i]; j++ ) { 46 printf(" %d\n",color->columns[i][j]); 47 } 48 printf("Number of rows %d\n",color->nrows[i]); 49 for ( j=0; j<color->nrows[i]; j++ ) { 50 printf(" %d %d \n",color->rows[i][j],color->columnsforrow[i][j]); 51 } 52 } 53 } 54 return 0; 55 } 56 57 /*@ 58 MatFDColoringSetParameters - Sets the parameters for the approximation of 59 Jacobian using finite differences. 60 61 $ J(u)_{:,i} = [J(u+h*dx_{i}) - J(u)]/h where 62 $ h = error_rel*u[i] if u[i] > umin 63 $ = error_rel*umin else 64 $ 65 $ dx_{i} = (0, ... 1, .... 0) 66 67 Input Parameters: 68 . coloring - the jacobian coloring context 69 . error_rel - relative error 70 . umin - minimum allowable u-value 71 72 .keywords: SNES, Jacobian, finite differences, parameters 73 @*/ 74 int MatFDColoringSetParameters(MatFDColoring matfd,double error,double umin) 75 { 76 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_COOKIE); 77 78 if (error != PETSC_DEFAULT) matfd->error_rel = error; 79 if (umin != PETSC_DEFAULT) matfd->umin = umin; 80 return 0; 81 } 82 83 /*@ 84 MatFDColoringSetFromOptions - Set coloring finite difference parameters from 85 the options database. 86 87 $ J(u)_{:,i} = [J(u+h*dx_{i}) - J(u)]/h where 88 $ h = error_rel*u[i] if u[i] > umin 89 $ = error_rel*.1 else 90 $ 91 $ dx_{i} = (0, ... 1, .... 0) 92 93 Input Parameters: 94 . coloring - the jacobian coloring context 95 96 Options Database: 97 . -mat_fd_coloring_error square root of relative error in function 98 . -mat_fd_coloring_umin see above 99 100 .keywords: SNES, Jacobian, finite differences, parameters 101 @*/ 102 int MatFDColoringSetFromOptions(MatFDColoring matfd) 103 { 104 int ierr,flag; 105 double error = PETSC_DEFAULT,umin = PETSC_DEFAULT; 106 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_COOKIE); 107 108 ierr = OptionsGetDouble(matfd->prefix,"-mat_fd_coloring_err",&error,&flag);CHKERRQ(ierr); 109 ierr = OptionsGetDouble(matfd->prefix,"-mat_fd_coloring_umin",&umin,&flag);CHKERRQ(ierr); 110 ierr = OptionsHasName(PETSC_NULL,"-help",&flag); CHKERRQ(ierr); 111 ierr = MatFDColoringSetParameters(matfd,error,umin); CHKERRQ(ierr); 112 if (flag) { 113 ierr = MatFDColoringPrintHelp(matfd); CHKERRQ(ierr); 114 } 115 return 0; 116 } 117 118 /*@ 119 MatFDColoringPrintHelp - Prints help message for matrix finite difference calculations 120 using coloring. 121 122 Input Parameter: 123 . fdcoloring - the MatFDColoring context 124 125 .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringSetFromOptions() 126 @*/ 127 int MatFDColoringPrintHelp(MatFDColoring fd) 128 { 129 PetscValidHeaderSpecific(fd,MAT_FDCOLORING_COOKIE); 130 131 PetscPrintf(fd->comm,"-mat_fd_coloring_err <err> set sqrt rel tol in function. Default 1.e-8\n"); 132 PetscPrintf(fd->comm,"-mat_fd_coloring_umin <umin> see users manual. Default 1.e-8\n"); 133 return 0; 134 } 135 136 /*@C 137 MatFDColoringCreate - Creates a matrix coloring context for finite difference 138 computation of Jacobians. 139 140 Input Parameters: 141 . mat - the matrix containing the nonzero structure of the Jacobian 142 . iscoloring - the coloring of the matrix 143 144 Output Parameter: 145 . color - the new coloring context 146 147 148 .seealso: MatFDColoringDestroy() 149 @*/ 150 int MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color) 151 { 152 MatFDColoring c; 153 MPI_Comm comm; 154 int ierr,M,N; 155 156 ierr = MatGetSize(mat,&M,&N); CHKERRQ(ierr); 157 if (M != N) SETERRQ(1,"MatFDColoringCreate:Only for square matrices"); 158 159 PetscObjectGetComm((PetscObject)mat,&comm); 160 PetscHeaderCreate(c,_MatFDColoring,MAT_FDCOLORING_COOKIE,0,comm); 161 PLogObjectCreate(c); 162 163 if (mat->ops.fdcoloringcreate) { 164 ierr = (*mat->ops.fdcoloringcreate)(mat,iscoloring,c); CHKERRQ(ierr); 165 } else { 166 SETERRQ(1,"MatFDColoringCreate:Code not yet written for this matrix type"); 167 } 168 169 c->error_rel = 1.e-8; 170 c->umin = 1.e-5; 171 172 *color = c; 173 174 return 0; 175 } 176 177 /*@C 178 MatFDColoringDestroy - Destroys a matrix coloring context that was created 179 via MatFDColoringCreate(). 180 181 Input Paramter: 182 . c - coloring context 183 184 .seealso: MatFDColoringCreate() 185 @*/ 186 int MatFDColoringDestroy(MatFDColoring c) 187 { 188 int i,ierr,flag; 189 190 ierr = OptionsHasName(PETSC_NULL,"-matfdcoloring_view",&flag); 191 if (flag) { 192 Viewer viewer; 193 ierr = ViewerFileOpenASCII(c->comm,"stdout",&viewer);CHKERRQ(ierr); 194 ierr = MatFDColoringView(c,viewer);CHKERRQ(ierr); 195 ierr = ViewerDestroy(viewer); CHKERRQ(ierr); 196 } 197 ierr = OptionsHasName(PETSC_NULL,"-matfdcoloring_view_info",&flag); 198 if (flag) { 199 Viewer viewer; 200 ierr = ViewerFileOpenASCII(c->comm,"stdout",&viewer);CHKERRQ(ierr); 201 ierr = ViewerPushFormat(viewer,VIEWER_FORMAT_ASCII_INFO, PETSC_NULL); CHKERRQ(ierr); 202 ierr = MatFDColoringView(c,viewer);CHKERRQ(ierr); 203 ierr = ViewerPopFormat(viewer); 204 ierr = ViewerDestroy(viewer); CHKERRQ(ierr); 205 } 206 207 for ( i=0; i<c->ncolors; i++ ) { 208 if (c->columns[i]) PetscFree(c->columns[i]); 209 if (c->rows[i]) PetscFree(c->rows[i]); 210 if (c->columnsforrow[i]) PetscFree(c->columnsforrow[i]); 211 } 212 PetscFree(c->ncolumns); 213 PetscFree(c->columns); 214 PetscFree(c->nrows); 215 PetscFree(c->rows); 216 PetscFree(c->columnsforrow); 217 PetscFree(c->scale); 218 PLogObjectDestroy(c); 219 PetscHeaderDestroy(c); 220 return 0; 221 } 222 223 /*@ 224 MatFDColoringApply - Given a matrix for which a MatFDColoring has been created, 225 computes the Jacobian for a function via finite differences. 226 227 Input Parameters: 228 . mat - location to store Jacobian 229 . coloring - coloring context created with MatFDColoringCreate() 230 . x1 - location at which Jacobian is to be computed 231 . w1,w2,w3 - three work vectors 232 . f - function for which Jacobian is required 233 . fctx - optional context required by function 234 235 .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView() 236 237 .keywords: coloring, Jacobian, finite differences 238 @*/ 239 int MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,Vec w1,Vec w2,Vec w3, 240 int (*f)(void *,Vec,Vec,void*),void *sctx,void *fctx) 241 { 242 int k, ierr,N,start,end,l,row,col,srow; 243 Scalar dx, mone = -1.0,*y,*scale = coloring->scale,*xx,*wscale = coloring->wscale; 244 double epsilon = coloring->error_rel, umin = coloring->umin; 245 MPI_Comm comm = coloring->comm; 246 247 ierr = MatZeroEntries(J); CHKERRQ(ierr); 248 249 ierr = VecGetOwnershipRange(x1,&start,&end); CHKERRQ(ierr); 250 ierr = VecGetSize(x1,&N); CHKERRQ(ierr); 251 ierr = VecGetArray(x1,&xx); CHKERRQ(ierr); 252 ierr = (*f)(sctx,x1,w1,fctx); CHKERRQ(ierr); 253 254 PetscMemzero(wscale,N*sizeof(Scalar)); 255 /* 256 Loop over each color 257 */ 258 259 for (k=0; k<coloring->ncolors; k++) { 260 ierr = VecCopy(x1,w3); CHKERRQ(ierr); 261 /* 262 Loop over each column associated with color adding the 263 perturbation to the vector w3. 264 */ 265 for (l=0; l<coloring->ncolumns[k]; l++) { 266 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 267 dx = xx[col-start]; 268 #if !defined(PETSC_COMPLEX) 269 if (dx < umin && dx >= 0.0) dx = .1; 270 else if (dx < 0.0 && dx > -umin) dx = -.1; 271 #else 272 if (abs(dx) < umin && real(dx) >= 0.0) dx = .1; 273 else if (real(dx) < 0.0 && abs(dx) < umin) dx = -.1; 274 #endif 275 dx *= epsilon; 276 wscale[col] = 1.0/dx; 277 VecSetValues(w3,1,&col,&dx,ADD_VALUES); 278 } 279 VecRestoreArray(x1,&xx); 280 /* 281 Evaluate function at x1 + dx (here dx is a vector, of perturbations) 282 */ 283 ierr = (*f)(sctx,w3,w2,fctx); CHKERRQ(ierr); 284 ierr = VecAXPY(&mone,w1,w2); CHKERRQ(ierr); 285 /* Communicate scale to all processors */ 286 #if !defined(PETSC_COMPLEX) 287 MPI_Allreduce(wscale,scale,N,MPI_DOUBLE,MPI_SUM,comm); 288 #else 289 MPI_Allreduce(wscale,scale,2*N,MPI_DOUBLE,MPI_SUM,comm); 290 #endif 291 /* 292 Loop over rows of vector putting results into Jacobian matrix 293 */ 294 VecGetArray(w2,&y); 295 for (l=0; l<coloring->nrows[k]; l++) { 296 row = coloring->rows[k][l]; 297 col = coloring->columnsforrow[k][l]; 298 y[row] *= scale[col]; 299 srow = row + start; 300 ierr = MatSetValues(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 301 } 302 VecRestoreArray(w2,&y); 303 } 304 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 305 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 306 return 0; 307 } 308