1 2 #ifndef lint 3 static char vcid[] = "$Id: fdmatrix.c,v 1.1 1996/10/09 17:52:05 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