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