Actual source code: fdmatrix.c
petsc-3.4.5 2014-06-29
2: /*
3: This is where the abstract matrix operations are defined that are
4: used for finite difference computations of Jacobians using coloring.
5: */
7: #include <petsc-private/matimpl.h> /*I "petscmat.h" I*/
11: PetscErrorCode MatFDColoringSetF(MatFDColoring fd,Vec F)
12: {
16: if (F) {
17: VecCopy(F,fd->w1);
18: fd->fset = PETSC_TRUE;
19: } else {
20: fd->fset = PETSC_FALSE;
21: }
22: return(0);
23: }
25: #include <petscdraw.h>
28: static PetscErrorCode MatFDColoringView_Draw_Zoom(PetscDraw draw,void *Aa)
29: {
30: MatFDColoring fd = (MatFDColoring)Aa;
32: PetscInt i,j;
33: PetscReal x,y;
36: /* loop over colors */
37: for (i=0; i<fd->ncolors; i++) {
38: for (j=0; j<fd->nrows[i]; j++) {
39: y = fd->M - fd->rows[i][j] - fd->rstart;
40: x = fd->columnsforrow[i][j];
41: PetscDrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1);
42: }
43: }
44: return(0);
45: }
49: static PetscErrorCode MatFDColoringView_Draw(MatFDColoring fd,PetscViewer viewer)
50: {
52: PetscBool isnull;
53: PetscDraw draw;
54: PetscReal xr,yr,xl,yl,h,w;
57: PetscViewerDrawGetDraw(viewer,0,&draw);
58: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
60: PetscObjectCompose((PetscObject)fd,"Zoomviewer",(PetscObject)viewer);
62: xr = fd->N; yr = fd->M; h = yr/10.0; w = xr/10.0;
63: xr += w; yr += h; xl = -w; yl = -h;
64: PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
65: PetscDrawZoom(draw,MatFDColoringView_Draw_Zoom,fd);
66: PetscObjectCompose((PetscObject)fd,"Zoomviewer",NULL);
67: return(0);
68: }
72: /*@C
73: MatFDColoringView - Views a finite difference coloring context.
75: Collective on MatFDColoring
77: Input Parameters:
78: + c - the coloring context
79: - viewer - visualization context
81: Level: intermediate
83: Notes:
84: The available visualization contexts include
85: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
86: . PETSC_VIEWER_STDOUT_WORLD - synchronized standard
87: output where only the first processor opens
88: the file. All other processors send their
89: data to the first processor to print.
90: - PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure
92: Notes:
93: Since PETSc uses only a small number of basic colors (currently 33), if the coloring
94: involves more than 33 then some seemingly identical colors are displayed making it look
95: like an illegal coloring. This is just a graphical artifact.
97: .seealso: MatFDColoringCreate()
99: .keywords: Mat, finite differences, coloring, view
100: @*/
101: PetscErrorCode MatFDColoringView(MatFDColoring c,PetscViewer viewer)
102: {
103: PetscErrorCode ierr;
104: PetscInt i,j;
105: PetscBool isdraw,iascii;
106: PetscViewerFormat format;
110: if (!viewer) {
111: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c),&viewer);
112: }
116: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
117: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
118: if (isdraw) {
119: MatFDColoringView_Draw(c,viewer);
120: } else if (iascii) {
121: PetscObjectPrintClassNamePrefixType((PetscObject)c,viewer,"MatFDColoring Object");
122: PetscViewerASCIIPrintf(viewer," Error tolerance=%G\n",c->error_rel);
123: PetscViewerASCIIPrintf(viewer," Umin=%G\n",c->umin);
124: PetscViewerASCIIPrintf(viewer," Number of colors=%D\n",c->ncolors);
126: PetscViewerGetFormat(viewer,&format);
127: if (format != PETSC_VIEWER_ASCII_INFO) {
128: for (i=0; i<c->ncolors; i++) {
129: PetscViewerASCIIPrintf(viewer," Information for color %D\n",i);
130: PetscViewerASCIIPrintf(viewer," Number of columns %D\n",c->ncolumns[i]);
131: for (j=0; j<c->ncolumns[i]; j++) {
132: PetscViewerASCIIPrintf(viewer," %D\n",c->columns[i][j]);
133: }
134: PetscViewerASCIIPrintf(viewer," Number of rows %D\n",c->nrows[i]);
135: for (j=0; j<c->nrows[i]; j++) {
136: PetscViewerASCIIPrintf(viewer," %D %D \n",c->rows[i][j],c->columnsforrow[i][j]);
137: }
138: }
139: }
140: PetscViewerFlush(viewer);
141: }
142: return(0);
143: }
147: /*@
148: MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
149: a Jacobian matrix using finite differences.
151: Logically Collective on MatFDColoring
153: The Jacobian is estimated with the differencing approximation
154: .vb
155: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
156: h = error_rel*u[i] if abs(u[i]) > umin
157: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
158: dx_{i} = (0, ... 1, .... 0)
159: .ve
161: Input Parameters:
162: + coloring - the coloring context
163: . error_rel - relative error
164: - umin - minimum allowable u-value magnitude
166: Level: advanced
168: .keywords: Mat, finite differences, coloring, set, parameters
170: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()
172: @*/
173: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
174: {
179: if (error != PETSC_DEFAULT) matfd->error_rel = error;
180: if (umin != PETSC_DEFAULT) matfd->umin = umin;
181: return(0);
182: }
188: /*@C
189: MatFDColoringGetFunction - Gets the function to use for computing the Jacobian.
191: Not Collective
193: Input Parameters:
194: . coloring - the coloring context
196: Output Parameters:
197: + f - the function
198: - fctx - the optional user-defined function context
200: Level: intermediate
202: .keywords: Mat, Jacobian, finite differences, set, function
204: .seealso: MatFDColoringCreate(), MatFDColoringSetFunction(), MatFDColoringSetFromOptions()
206: @*/
207: PetscErrorCode MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
208: {
211: if (f) *f = matfd->f;
212: if (fctx) *fctx = matfd->fctx;
213: return(0);
214: }
218: /*@C
219: MatFDColoringSetFunction - Sets the function to use for computing the Jacobian.
221: Logically Collective on MatFDColoring
223: Input Parameters:
224: + coloring - the coloring context
225: . f - the function
226: - fctx - the optional user-defined function context
228: Calling sequence of (*f) function:
229: For SNES: PetscErrorCode (*f)(SNES,Vec,Vec,void*)
230: If not using SNES: PetscErrorCode (*f)(void *dummy,Vec,Vec,void*) and dummy is ignored
232: Level: advanced
234: Notes: This function is usually used automatically by SNES (when one uses SNESSetJacobian() with the argument
235: SNESComputeJacobianDefaultColor()) and only needs to be used by someone computing a matrix via coloring directly by
236: calling MatFDColoringApply()
238: Fortran Notes:
239: In Fortran you must call MatFDColoringSetFunction() for a coloring object to
240: be used without SNES or within the SNES solvers.
242: .keywords: Mat, Jacobian, finite differences, set, function
244: .seealso: MatFDColoringCreate(), MatFDColoringGetFunction(), MatFDColoringSetFromOptions()
246: @*/
247: PetscErrorCode MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
248: {
251: matfd->f = f;
252: matfd->fctx = fctx;
253: return(0);
254: }
258: /*@
259: MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
260: the options database.
262: Collective on MatFDColoring
264: The Jacobian, F'(u), is estimated with the differencing approximation
265: .vb
266: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
267: h = error_rel*u[i] if abs(u[i]) > umin
268: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
269: dx_{i} = (0, ... 1, .... 0)
270: .ve
272: Input Parameter:
273: . coloring - the coloring context
275: Options Database Keys:
276: + -mat_fd_coloring_err <err> - Sets <err> (square root
277: of relative error in the function)
278: . -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
279: . -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
280: . -mat_fd_coloring_view - Activates basic viewing
281: . -mat_fd_coloring_view ::ascii_info - Activates viewing info
282: - -mat_fd_coloring_view draw - Activates drawing
284: Level: intermediate
286: .keywords: Mat, finite differences, parameters
288: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()
290: @*/
291: PetscErrorCode MatFDColoringSetFromOptions(MatFDColoring matfd)
292: {
294: PetscBool flg;
295: char value[3];
300: PetscObjectOptionsBegin((PetscObject)matfd);
301: PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,0);
302: PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,0);
303: PetscOptionsString("-mat_fd_type","Algorithm to compute h, wp or ds","MatFDColoringCreate",matfd->htype,value,3,&flg);
304: if (flg) {
305: if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
306: else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
307: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",value);
308: }
309: /* process any options handlers added with PetscObjectAddOptionsHandler() */
310: PetscObjectProcessOptionsHandlers((PetscObject)matfd);
311: PetscOptionsEnd();
312: return(0);
313: }
317: PetscErrorCode MatFDColoringViewFromOptions(MatFDColoring fd,const char prefix[],const char optionname[])
318: {
319: PetscErrorCode ierr;
320: PetscBool flg;
321: PetscViewer viewer;
322: PetscViewerFormat format;
325: if (prefix) {
326: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),prefix,optionname,&viewer,&format,&flg);
327: } else {
328: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->prefix,optionname,&viewer,&format,&flg);
329: }
330: if (flg) {
331: PetscViewerPushFormat(viewer,format);
332: MatFDColoringView(fd,viewer);
333: PetscViewerPopFormat(viewer);
334: PetscViewerDestroy(&viewer);
335: }
336: return(0);
337: }
341: /*@
342: MatFDColoringCreate - Creates a matrix coloring context for finite difference
343: computation of Jacobians.
345: Collective on Mat
347: Input Parameters:
348: + mat - the matrix containing the nonzero structure of the Jacobian
349: - iscoloring - the coloring of the matrix; usually obtained with MatGetColoring() or DMCreateColoring()
351: Output Parameter:
352: . color - the new coloring context
354: Level: intermediate
356: .seealso: MatFDColoringDestroy(),SNESComputeJacobianDefaultColor(), ISColoringCreate(),
357: MatFDColoringSetFunction(), MatFDColoringSetFromOptions(), MatFDColoringApply(),
358: MatFDColoringView(), MatFDColoringSetParameters(), MatGetColoring(), DMCreateColoring()
359: @*/
360: PetscErrorCode MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color)
361: {
362: MatFDColoring c;
363: MPI_Comm comm;
365: PetscInt M,N;
368: PetscLogEventBegin(MAT_FDColoringCreate,mat,0,0,0);
369: MatGetSize(mat,&M,&N);
370: if (M != N) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Only for square matrices");
371: PetscObjectGetComm((PetscObject)mat,&comm);
372: PetscHeaderCreate(c,_p_MatFDColoring,int,MAT_FDCOLORING_CLASSID,"MatFDColoring","Jacobian computation via finite differences with coloring","Mat",comm,MatFDColoringDestroy,MatFDColoringView);
374: c->ctype = iscoloring->ctype;
376: if (mat->ops->fdcoloringcreate) {
377: (*mat->ops->fdcoloringcreate)(mat,iscoloring,c);
378: } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);
380: MatGetVecs(mat,NULL,&c->w1);
381: PetscLogObjectParent(c,c->w1);
382: VecDuplicate(c->w1,&c->w2);
383: PetscLogObjectParent(c,c->w2);
385: c->error_rel = PETSC_SQRT_MACHINE_EPSILON;
386: c->umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
387: c->currentcolor = -1;
388: c->htype = "wp";
389: c->fset = PETSC_FALSE;
391: *color = c;
392: PetscObjectCompose((PetscObject)mat,"SNESMatFDColoring",(PetscObject)c);
393: PetscLogEventEnd(MAT_FDColoringCreate,mat,0,0,0);
394: return(0);
395: }
399: /*@
400: MatFDColoringDestroy - Destroys a matrix coloring context that was created
401: via MatFDColoringCreate().
403: Collective on MatFDColoring
405: Input Parameter:
406: . c - coloring context
408: Level: intermediate
410: .seealso: MatFDColoringCreate()
411: @*/
412: PetscErrorCode MatFDColoringDestroy(MatFDColoring *c)
413: {
415: PetscInt i;
418: if (!*c) return(0);
419: if (--((PetscObject)(*c))->refct > 0) {*c = 0; return(0);}
421: for (i=0; i<(*c)->ncolors; i++) {
422: PetscFree((*c)->columns[i]);
423: PetscFree((*c)->rows[i]);
424: PetscFree((*c)->columnsforrow[i]);
425: if ((*c)->vscaleforrow) {PetscFree((*c)->vscaleforrow[i]);}
426: }
427: PetscFree((*c)->ncolumns);
428: PetscFree((*c)->columns);
429: PetscFree((*c)->nrows);
430: PetscFree((*c)->rows);
431: PetscFree((*c)->columnsforrow);
432: PetscFree((*c)->vscaleforrow);
433: VecDestroy(&(*c)->vscale);
434: VecDestroy(&(*c)->w1);
435: VecDestroy(&(*c)->w2);
436: VecDestroy(&(*c)->w3);
437: PetscHeaderDestroy(c);
438: return(0);
439: }
443: /*@C
444: MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
445: that are currently being perturbed.
447: Not Collective
449: Input Parameters:
450: . coloring - coloring context created with MatFDColoringCreate()
452: Output Parameters:
453: + n - the number of local columns being perturbed
454: - cols - the column indices, in global numbering
456: Level: intermediate
458: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringApply()
460: .keywords: coloring, Jacobian, finite differences
461: @*/
462: PetscErrorCode MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt *n,PetscInt *cols[])
463: {
465: if (coloring->currentcolor >= 0) {
466: *n = coloring->ncolumns[coloring->currentcolor];
467: *cols = coloring->columns[coloring->currentcolor];
468: } else {
469: *n = 0;
470: }
471: return(0);
472: }
476: /*@
477: MatFDColoringApply - Given a matrix for which a MatFDColoring context
478: has been created, computes the Jacobian for a function via finite differences.
480: Collective on MatFDColoring
482: Input Parameters:
483: + mat - location to store Jacobian
484: . coloring - coloring context created with MatFDColoringCreate()
485: . x1 - location at which Jacobian is to be computed
486: - sctx - context required by function, if this is being used with the SNES solver then it is SNES object, otherwise it is null
488: Options Database Keys:
489: + -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
490: . -mat_fd_coloring_view - Activates basic viewing or coloring
491: . -mat_fd_coloring_view draw - Activates drawing of coloring
492: - -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info
494: Level: intermediate
496: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringSetFunction()
498: .keywords: coloring, Jacobian, finite differences
499: @*/
500: PetscErrorCode MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx)
501: {
508: if (!coloring->f) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");
509: if (!J->ops->fdcoloringapply) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not supported for this matrix type %s",((PetscObject)J)->type_name);
510: (*J->ops->fdcoloringapply)(J,coloring,x1,flag,sctx);
511: return(0);
512: }
516: PetscErrorCode MatFDColoringApply_AIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx)
517: {
518: PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void*))coloring->f;
520: PetscInt k,start,end,l,row,col,srow,**vscaleforrow;
521: PetscScalar dx,*y,*xx,*w3_array;
522: PetscScalar *vscale_array;
523: PetscReal epsilon = coloring->error_rel,umin = coloring->umin,unorm;
524: Vec w1 = coloring->w1,w2=coloring->w2,w3;
525: void *fctx = coloring->fctx;
526: PetscBool flg = PETSC_FALSE;
527: PetscInt ctype = coloring->ctype,N,col_start=0,col_end=0;
528: Vec x1_tmp;
534: if (!f) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");
536: PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);
537: MatSetUnfactored(J);
538: PetscOptionsGetBool(NULL,"-mat_fd_coloring_dont_rezero",&flg,NULL);
539: if (flg) {
540: PetscInfo(coloring,"Not calling MatZeroEntries()\n");
541: } else {
542: PetscBool assembled;
543: MatAssembled(J,&assembled);
544: if (assembled) {
545: MatZeroEntries(J);
546: }
547: }
549: x1_tmp = x1;
550: if (!coloring->vscale) {
551: VecDuplicate(x1_tmp,&coloring->vscale);
552: }
554: if (coloring->htype[0] == 'w') { /* tacky test; need to make systematic if we add other approaches to computing h*/
555: VecNorm(x1_tmp,NORM_2,&unorm);
556: }
557: VecGetOwnershipRange(w1,&start,&end); /* OwnershipRange is used by ghosted x! */
559: /* Set w1 = F(x1) */
560: if (!coloring->fset) {
561: PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
562: (*f)(sctx,x1_tmp,w1,fctx);
563: PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
564: } else {
565: coloring->fset = PETSC_FALSE;
566: }
568: if (!coloring->w3) {
569: VecDuplicate(x1_tmp,&coloring->w3);
570: PetscLogObjectParent(coloring,coloring->w3);
571: }
572: w3 = coloring->w3;
574: /* Compute all the local scale factors, including ghost points */
575: VecGetLocalSize(x1_tmp,&N);
576: VecGetArray(x1_tmp,&xx);
577: VecGetArray(coloring->vscale,&vscale_array);
578: if (ctype == IS_COLORING_GHOSTED) {
579: col_start = 0; col_end = N;
580: } else if (ctype == IS_COLORING_GLOBAL) {
581: xx = xx - start;
582: vscale_array = vscale_array - start;
583: col_start = start; col_end = N + start;
584: }
585: for (col=col_start; col<col_end; col++) {
586: /* Loop over each local column, vscale[col] = 1./(epsilon*dx[col]) */
587: if (coloring->htype[0] == 'w') {
588: dx = 1.0 + unorm;
589: } else {
590: dx = xx[col];
591: }
592: if (dx == (PetscScalar)0.0) dx = 1.0;
593: if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin;
594: else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
595: dx *= epsilon;
596: vscale_array[col] = (PetscScalar)1.0/dx;
597: }
598: if (ctype == IS_COLORING_GLOBAL) vscale_array = vscale_array + start;
599: VecRestoreArray(coloring->vscale,&vscale_array);
600: if (ctype == IS_COLORING_GLOBAL) {
601: VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
602: VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);
603: }
605: if (coloring->vscaleforrow) {
606: vscaleforrow = coloring->vscaleforrow;
607: } else SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_NULL,"Null Object: coloring->vscaleforrow");
609: /*
610: Loop over each color
611: */
612: VecGetArray(coloring->vscale,&vscale_array);
613: for (k=0; k<coloring->ncolors; k++) {
614: coloring->currentcolor = k;
616: VecCopy(x1_tmp,w3);
617: VecGetArray(w3,&w3_array);
618: if (ctype == IS_COLORING_GLOBAL) w3_array = w3_array - start;
619: /*
620: Loop over each column associated with color
621: adding the perturbation to the vector w3.
622: */
623: for (l=0; l<coloring->ncolumns[k]; l++) {
624: col = coloring->columns[k][l]; /* local column of the matrix we are probing for */
625: if (coloring->htype[0] == 'w') {
626: dx = 1.0 + unorm;
627: } else {
628: dx = xx[col];
629: }
630: if (dx == (PetscScalar)0.0) dx = 1.0;
631: if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin;
632: else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
633: dx *= epsilon;
634: if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter");
635: w3_array[col] += dx;
636: }
637: if (ctype == IS_COLORING_GLOBAL) w3_array = w3_array + start;
638: VecRestoreArray(w3,&w3_array);
640: /*
641: Evaluate function at w3 = x1 + dx (here dx is a vector of perturbations)
642: w2 = F(x1 + dx) - F(x1)
643: */
644: PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
645: (*f)(sctx,w3,w2,fctx);
646: PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
647: VecAXPY(w2,-1.0,w1);
649: /*
650: Loop over rows of vector, putting results into Jacobian matrix
651: */
652: VecGetArray(w2,&y);
653: for (l=0; l<coloring->nrows[k]; l++) {
654: row = coloring->rows[k][l]; /* local row index */
655: col = coloring->columnsforrow[k][l]; /* global column index */
656: y[row] *= vscale_array[vscaleforrow[k][l]];
657: srow = row + start;
658: MatSetValues(J,1,&srow,1,&col,y+row,INSERT_VALUES);
659: }
660: VecRestoreArray(w2,&y);
661: } /* endof for each color */
662: if (ctype == IS_COLORING_GLOBAL) xx = xx + start;
663: VecRestoreArray(coloring->vscale,&vscale_array);
664: VecRestoreArray(x1_tmp,&xx);
666: coloring->currentcolor = -1;
668: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
669: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
670: PetscLogEventEnd(MAT_FDColoringApply,coloring,J,x1,0);
672: MatFDColoringViewFromOptions(coloring,NULL,"-mat_fd_coloring_view");
673: return(0);
674: }