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