1 2 #ifndef lint 3 static char vcid[] = "$Id: snesj2.c,v 1.3 1996/08/27 18:38:51 bsmith Exp bsmith $"; 4 #endif 5 6 #include "snesimpl.h" /*I "snes.h" I*/ 7 8 #include "src/mat/impls/aij/seq/aij.h" 9 10 11 /* ---------------------------------------------------------------------------------*/ 12 13 14 /*@C 15 SNESDefaultComputeJacobianWithColoring 16 17 Input Parameters: 18 . snes - nonlinear solver object 19 . x1 - location at which to evaluate Jacobian 20 . ctx - MatFDColoring contex 21 22 Output Parameters: 23 . J - Jacobian matrix 24 . B - Jacobian preconditioner 25 . flag - flag indicating if the matrix nonzero structure has changed 26 27 .keywords: SNES, finite differences, Jacobian 28 29 .seealso: SNESSetJacobian(), SNESTestJacobian() 30 @*/ 31 int SNESDefaultComputeJacobianWithColoring(SNES snes,Vec x1,Mat *JJ,Mat *B,MatStructure *flag,void *ctx) 32 { 33 MatFDColoring color = (MatFDColoring) ctx; 34 Mat J = *JJ; 35 Vec jj1,jj2,x2; 36 int k, ierr,N,start,end,l,row,col; 37 Scalar dx, mone = -1.0,*y,*scale = color->scale,*xx,*wscale = color->wscale; 38 double epsilon = 1.e-8; /* assumes double precision */ 39 MPI_Comm comm = color->comm; 40 41 PetscTrValid(0,0); 42 ierr = MatZeroEntries(J); CHKERRQ(ierr); 43 if (!snes->nvwork) { 44 ierr = VecDuplicateVecs(x1,3,&snes->vwork); CHKERRQ(ierr); 45 snes->nvwork = 3; 46 PLogObjectParents(snes,3,snes->vwork); 47 } 48 jj1 = snes->vwork[0]; jj2 = snes->vwork[1]; x2 = snes->vwork[2]; 49 50 ierr = VecGetOwnershipRange(x1,&start,&end); CHKERRQ(ierr); 51 ierr = VecGetSize(x1,&N); CHKERRQ(ierr); 52 ierr = VecGetArray(x1,&xx); CHKERRQ(ierr); 53 ierr = SNESComputeFunction(snes,x1,jj1); CHKERRQ(ierr); 54 55 PetscMemzero(wscale,N*sizeof(Scalar)); 56 /* 57 Loop over each color 58 */ 59 for (k=0; k<color->ncolors; k++) { 60 ierr = VecCopy(x1,x2); CHKERRQ(ierr); 61 /* 62 Loop over each column associated with color adding the 63 perturbation to the vector x2. 64 */ 65 for (l=0; l<color->ncolumns[k]; l++) { 66 col = color->columns[k][l]; /* column of the matrix we are probing for */ 67 dx = xx[col-start]; 68 #if !defined(PETSC_COMPLEX) 69 if (dx < 1.e-16 && dx >= 0.0) dx = 1.e-1; 70 else if (dx < 0.0 && dx > -1.e-16) dx = -1.e-1; 71 #else 72 if (abs(dx) < 1.e-16 && real(dx) >= 0.0) dx = 1.e-1; 73 else if (real(dx) < 0.0 && abs(dx) < 1.e-16) dx = -1.e-1; 74 #endif 75 dx *= epsilon; 76 wscale[col] = 1.0/dx; 77 VecSetValues(x2,1,&col,&dx,ADD_VALUES); 78 } 79 VecRestoreArray(x1,&xx); 80 /* 81 Evaluate function at x1 + dx (here dx is a vector, of perturbations) 82 */ 83 ierr = SNESComputeFunction(snes,x2,jj2); CHKERRQ(ierr); 84 ierr = VecAXPY(&mone,jj1,jj2); CHKERRQ(ierr); 85 86 /* Communicate scale to all processors */ 87 #if !defined(PETSC_COMPLEX) 88 MPI_Allreduce(wscale,scale,N,MPI_DOUBLE,MPI_SUM,comm); 89 #else 90 MPI_Allreduce(wscale,scale,2*N,MPI_DOUBLE,MPI_SUM,comm); 91 #endif 92 93 /* 94 Loop over rows of vector putting results into Jacobian matrix 95 */ 96 VecGetArray(jj2,&y); 97 for (l=0; l<color->nrows[k]; l++) { 98 row = color->rows[k][l]; 99 col = color->columnsforrow[k][l]; 100 y[row-start] *= scale[col]; 101 ierr = MatSetValues(J,1,&row,1,&col,y+row-start,INSERT_VALUES);CHKERRQ(ierr); 102 } 103 VecRestoreArray(jj2,&y); 104 } 105 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 106 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 107 *flag = SAME_NONZERO_PATTERN; 108 return 0; 109 } 110 111 112 113