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