Actual source code: pack.c
1: #define PETSCDM_DLL
2:
3: #include petscda.h
4: #include src/dm/da/daimpl.h
5: #include petscmat.h
7: typedef struct _VecPackOps *VecPackOps;
8: struct _VecPackOps {
9: PetscErrorCode (*view)(VecPack,PetscViewer);
10: PetscErrorCode (*createglobalvector)(VecPack,Vec*);
11: PetscErrorCode (*getcoloring)(VecPack,ISColoringType,ISColoring*);
12: PetscErrorCode (*getmatrix)(VecPack, MatType,Mat*);
13: PetscErrorCode (*getinterpolation)(VecPack,VecPack,Mat*,Vec*);
14: PetscErrorCode (*refine)(VecPack,MPI_Comm,VecPack*);
15: PetscErrorCode (*getinjection)(VecPack,VecPack,VecScatter*);
17: PetscErrorCode (*forminitialguess)(VecPack,PetscErrorCode (*)(void),Vec,void*);
18: PetscErrorCode (*formfunction)(VecPack,PetscErrorCode (*)(void),Vec,void*);
19: };
21: /*
22: rstart is where an array/subvector starts in the global parallel vector, so arrays
23: rstarts are meaningless (and set to the previous one) except on processor 0
24: */
26: typedef enum {VECPACK_ARRAY, VECPACK_DA, VECPACK_VECSCATTER} VecPackLinkType;
28: struct VecPackLink {
29: VecPackLinkType type;
30: struct VecPackLink *next;
31: PetscInt n,rstart; /* rstart is relative to this processor */
33: /* only used for VECPACK_DA */
34: PetscInt *grstarts; /* global row for first unknown of this DA on each process */
35: DA da;
37: /* only used for VECPACK_ARRAY */
38: PetscInt grstart; /* global row for first array unknown */
39: PetscMPIInt rank; /* process where array unknowns live */
40: };
42: struct _p_VecPack {
43: PETSCHEADER(struct _VecPackOps);
44: PetscInt n,N,rstart; /* rstart is relative to all processors, n unknowns owned by this process, N is total unknowns */
45: PetscInt nghost; /* number of all local entries include DA ghost points and any shared redundant arrays */
46: Vec globalvector;
47: PetscInt nDA,nredundant; /* how many DA's and seperate redundant arrays used to build VecPack */
48: struct VecPackLink *next;
49: };
53: /*
54: Maps from
55: */
56: PetscErrorCode VecPackFormInitialGuess_DADADA(VecPack pack,PetscErrorCode (*fun)(void),Vec X,void *ctx)
57: {
59: PetscErrorCode (*f)(DALocalInfo*,void*,DALocalInfo*,void*,DALocalInfo*,void*) =
60: (PetscErrorCode (*)(DALocalInfo*,void*,DALocalInfo*,void*,DALocalInfo*,void*)) fun;
61: DALocalInfo da1,da2,da3;
62: DA DA1,DA2,DA3;
63: void *x1,*x2,*x3;
64: Vec X1,X2,X3;
67: VecPackGetEntries(pack,&DA1,&DA2,&DA3);
68: DAGetLocalInfo(DA1,&da1);
69: DAGetLocalInfo(DA2,&da2);
70: DAGetLocalInfo(DA3,&da3);
71: VecPackGetAccess(pack,X,&X1,&X2,&X3);
72: DAVecGetArray(DA1,X1,&x1);
73: DAVecGetArray(DA2,X2,&x2);
74: DAVecGetArray(DA3,X3,&x3);
77: (*f)(&da1,x1,&da2,x2,&da3,x3);
79: return(0);
80: }
84: PetscErrorCode VecView_DADADA(VecPack pack,Vec X,PetscViewer viewer)
85: {
87: DA DA1,DA2,DA3;
88: Vec X1,X2,X3;
91: VecPackGetEntries(pack,&DA1,&DA2,&DA3);
92: VecPackGetAccess(pack,X,&X1,&X2,&X3);
93: VecView(X1,viewer);
94: VecView(X2,viewer);
95: VecView(X3,viewer);
96: return(0);
97: }
101: /*@C
102: VecPackCreate - Creates a vector packer, used to generate "composite"
103: vectors made up of several subvectors.
105: Collective on MPI_Comm
107: Input Parameter:
108: . comm - the processors that will share the global vector
110: Output Parameters:
111: . packer - the packer object
113: Level: advanced
115: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackScatter(),
116: VecPackGather(), VecPackCreateGlobalVector(), VecPackGetGlobalIndices(), VecPackGetAccess()
117: VecPackGetLocalVectors(), VecPackRestoreLocalVectors()
119: @*/
120: PetscErrorCode VecPackCreate(MPI_Comm comm,VecPack *packer)
121: {
123: VecPack p;
127: *packer = PETSC_NULL;
128: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
129: DMInitializePackage(PETSC_NULL);
130: #endif
132: PetscHeaderCreate(p,_p_VecPack,struct _VecPackOps,DA_COOKIE,0,"VecPack",comm,VecPackDestroy,0);
133: p->n = 0;
134: p->next = PETSC_NULL;
135: p->comm = comm;
136: p->globalvector = PETSC_NULL;
137: p->nredundant = 0;
138: p->nDA = 0;
140: p->ops->createglobalvector = VecPackCreateGlobalVector;
141: p->ops->refine = VecPackRefine;
142: p->ops->getinterpolation = VecPackGetInterpolation;
143: p->ops->getmatrix = VecPackGetMatrix;
144: p->ops->getcoloring = VecPackGetColoring;
146: p->ops->forminitialguess = VecPackFormInitialGuess_DADADA;
147: *packer = p;
148: return(0);
149: }
153: /*@C
154: VecPackDestroy - Destroys a vector packer.
156: Collective on VecPack
158: Input Parameter:
159: . packer - the packer object
161: Level: advanced
163: .seealso VecPackCreate(), VecPackAddArray(), VecPackAddDA(), VecPackScatter(),
164: VecPackGather(), VecPackCreateGlobalVector(), VecPackGetGlobalIndices(), VecPackGetAccess()
166: @*/
167: PetscErrorCode VecPackDestroy(VecPack packer)
168: {
169: PetscErrorCode ierr;
170: struct VecPackLink *next = packer->next,*prev;
173: if (--packer->refct > 0) return(0);
174: while (next) {
175: prev = next;
176: next = next->next;
177: if (prev->type == VECPACK_DA) {
178: DADestroy(prev->da);
179: }
180: if (prev->grstarts) {
181: PetscFree(prev->grstarts);
182: }
183: PetscFree(prev);
184: }
185: if (packer->globalvector) {
186: VecDestroy(packer->globalvector);
187: }
188: PetscHeaderDestroy(packer);
189: return(0);
190: }
192: /* --------------------------------------------------------------------------------------*/
196: PetscErrorCode VecPackGetAccess_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar **array)
197: {
199: PetscScalar *varray;
200: PetscMPIInt rank;
203: MPI_Comm_rank(packer->comm,&rank);
204: if (array) {
205: if (rank == mine->rank) {
206: VecGetArray(vec,&varray);
207: *array = varray + mine->rstart;
208: VecRestoreArray(vec,&varray);
209: } else {
210: *array = 0;
211: }
212: }
213: return(0);
214: }
218: PetscErrorCode VecPackGetAccess_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec *global)
219: {
221: PetscScalar *array;
224: if (global) {
225: DAGetGlobalVector(mine->da,global);
226: VecGetArray(vec,&array);
227: VecPlaceArray(*global,array+mine->rstart);
228: VecRestoreArray(vec,&array);
229: }
230: return(0);
231: }
235: PetscErrorCode VecPackRestoreAccess_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar **array)
236: {
238: return(0);
239: }
243: PetscErrorCode VecPackRestoreAccess_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec *global)
244: {
248: if (global) {
249: VecResetArray(*global);
250: DARestoreGlobalVector(mine->da,global);
251: }
252: return(0);
253: }
257: PetscErrorCode VecPackScatter_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar *array)
258: {
260: PetscScalar *varray;
261: PetscMPIInt rank;
264: MPI_Comm_rank(packer->comm,&rank);
265: if (rank == mine->rank) {
266: VecGetArray(vec,&varray);
267: PetscMemcpy(array,varray+mine->rstart,mine->n*sizeof(PetscScalar));
268: VecRestoreArray(vec,&varray);
269: }
270: MPI_Bcast(array,mine->n,MPIU_SCALAR,mine->rank,packer->comm);
271: return(0);
272: }
276: PetscErrorCode VecPackScatter_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec local)
277: {
279: PetscScalar *array;
280: Vec global;
283: DAGetGlobalVector(mine->da,&global);
284: VecGetArray(vec,&array);
285: VecPlaceArray(global,array+mine->rstart);
286: DAGlobalToLocalBegin(mine->da,global,INSERT_VALUES,local);
287: DAGlobalToLocalEnd(mine->da,global,INSERT_VALUES,local);
288: VecRestoreArray(vec,&array);
289: VecResetArray(global);
290: DARestoreGlobalVector(mine->da,&global);
291: return(0);
292: }
296: PetscErrorCode VecPackGather_Array(VecPack packer,struct VecPackLink *mine,Vec vec,PetscScalar *array)
297: {
299: PetscScalar *varray;
300: PetscMPIInt rank;
303: MPI_Comm_rank(packer->comm,&rank);
304: if (rank == mine->rank) {
305: VecGetArray(vec,&varray);
306: if (varray+mine->rstart == array) SETERRQ(PETSC_ERR_ARG_WRONG,"You need not VecPackGather() into objects obtained via VecPackGetAccess()");
307: PetscMemcpy(varray+mine->rstart,array,mine->n*sizeof(PetscScalar));
308: VecRestoreArray(vec,&varray);
309: }
310: return(0);
311: }
315: PetscErrorCode VecPackGather_DA(VecPack packer,struct VecPackLink *mine,Vec vec,Vec local)
316: {
318: PetscScalar *array;
319: Vec global;
322: DAGetGlobalVector(mine->da,&global);
323: VecGetArray(vec,&array);
324: VecPlaceArray(global,array+mine->rstart);
325: DALocalToGlobal(mine->da,local,INSERT_VALUES,global);
326: VecRestoreArray(vec,&array);
327: VecResetArray(global);
328: DARestoreGlobalVector(mine->da,&global);
329: return(0);
330: }
332: /* ----------------------------------------------------------------------------------*/
334: #include <stdarg.h>
338: /*@C
339: VecPackGetAccess - Allows one to access the individual packed vectors in their global
340: representation.
342: Collective on VecPack
344: Input Parameter:
345: + packer - the packer object
346: . gvec - the global vector
347: - ... - the individual sequential or parallel objects (arrays or vectors)
348:
349: Level: advanced
351: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
352: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackScatter(),
353: VecPackRestoreAccess(), VecPackGetLocalVectors(), VecPackRestoreLocalVectors()
355: @*/
356: PetscErrorCode VecPackGetAccess(VecPack packer,Vec gvec,...)
357: {
358: va_list Argp;
359: PetscErrorCode ierr;
360: struct VecPackLink *next = packer->next;
363: if (!packer->globalvector) {
364: SETERRQ(PETSC_ERR_ORDER,"Must first create global vector with VecPackCreateGlobalVector()");
365: }
367: /* loop over packed objects, handling one at at time */
368: va_start(Argp,gvec);
369: while (next) {
370: if (next->type == VECPACK_ARRAY) {
371: PetscScalar **array;
372: array = va_arg(Argp, PetscScalar**);
373: VecPackGetAccess_Array(packer,next,gvec,array);
374: } else if (next->type == VECPACK_DA) {
375: Vec *vec;
376: vec = va_arg(Argp, Vec*);
377: VecPackGetAccess_DA(packer,next,gvec,vec);
378: } else {
379: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
380: }
381: next = next->next;
382: }
383: va_end(Argp);
384: return(0);
385: }
389: /*@C
390: VecPackRestoreAccess - Allows one to access the individual packed vectors in their global
391: representation.
393: Collective on VecPack
395: Input Parameter:
396: + packer - the packer object
397: . gvec - the global vector
398: - ... - the individual sequential or parallel objects (arrays or vectors)
399:
400: Level: advanced
402: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
403: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackScatter(),
404: VecPackRestoreAccess()
406: @*/
407: PetscErrorCode VecPackRestoreAccess(VecPack packer,Vec gvec,...)
408: {
409: va_list Argp;
410: PetscErrorCode ierr;
411: struct VecPackLink *next = packer->next;
414: if (!packer->globalvector) {
415: SETERRQ(PETSC_ERR_ORDER,"Must first create global vector with VecPackCreateGlobalVector()");
416: }
418: /* loop over packed objects, handling one at at time */
419: va_start(Argp,gvec);
420: while (next) {
421: if (next->type == VECPACK_ARRAY) {
422: PetscScalar **array;
423: array = va_arg(Argp, PetscScalar**);
424: VecPackRestoreAccess_Array(packer,next,gvec,array);
425: } else if (next->type == VECPACK_DA) {
426: Vec *vec;
427: vec = va_arg(Argp, Vec*);
428: VecPackRestoreAccess_DA(packer,next,gvec,vec);
429: } else {
430: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
431: }
432: next = next->next;
433: }
434: va_end(Argp);
435: return(0);
436: }
440: /*@C
441: VecPackScatter - Scatters from a global packed vector into its individual local vectors
443: Collective on VecPack
445: Input Parameter:
446: + packer - the packer object
447: . gvec - the global vector
448: - ... - the individual sequential objects (arrays or vectors)
449:
450: Level: advanced
452: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
453: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
454: VecPackGetLocalVectors(), VecPackRestoreLocalVectors()
456: @*/
457: PetscErrorCode VecPackScatter(VecPack packer,Vec gvec,...)
458: {
459: va_list Argp;
460: PetscErrorCode ierr;
461: struct VecPackLink *next = packer->next;
464: if (!packer->globalvector) {
465: SETERRQ(PETSC_ERR_ORDER,"Must first create global vector with VecPackCreateGlobalVector()");
466: }
468: /* loop over packed objects, handling one at at time */
469: va_start(Argp,gvec);
470: while (next) {
471: if (next->type == VECPACK_ARRAY) {
472: PetscScalar *array;
473: array = va_arg(Argp, PetscScalar*);
474: VecPackScatter_Array(packer,next,gvec,array);
475: } else if (next->type == VECPACK_DA) {
476: Vec vec;
477: vec = va_arg(Argp, Vec);
479: VecPackScatter_DA(packer,next,gvec,vec);
480: } else {
481: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
482: }
483: next = next->next;
484: }
485: va_end(Argp);
486: return(0);
487: }
491: /*@C
492: VecPackGather - Gathers into a global packed vector from its individual local vectors
494: Collective on VecPack
496: Input Parameter:
497: + packer - the packer object
498: . gvec - the global vector
499: - ... - the individual sequential objects (arrays or vectors)
500:
501: Level: advanced
503: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
504: VecPackScatter(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
505: VecPackGetLocalVectors(), VecPackRestoreLocalVectors()
507: @*/
508: PetscErrorCode VecPackGather(VecPack packer,Vec gvec,...)
509: {
510: va_list Argp;
511: PetscErrorCode ierr;
512: struct VecPackLink *next = packer->next;
515: if (!packer->globalvector) {
516: SETERRQ(PETSC_ERR_ORDER,"Must first create global vector with VecPackCreateGlobalVector()");
517: }
519: /* loop over packed objects, handling one at at time */
520: va_start(Argp,gvec);
521: while (next) {
522: if (next->type == VECPACK_ARRAY) {
523: PetscScalar *array;
524: array = va_arg(Argp, PetscScalar*);
525: VecPackGather_Array(packer,next,gvec,array);
526: } else if (next->type == VECPACK_DA) {
527: Vec vec;
528: vec = va_arg(Argp, Vec);
530: VecPackGather_DA(packer,next,gvec,vec);
531: } else {
532: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
533: }
534: next = next->next;
535: }
536: va_end(Argp);
537: return(0);
538: }
542: /*@C
543: VecPackAddArray - adds an "redundant" array to a VecPack. The array values will
544: be stored in part of the array on process orank.
546: Collective on VecPack
548: Input Parameter:
549: + packer - the packer object
550: . orank - the process on which the array entries officially live, this number must be
551: the same on all processes.
552: - n - the length of the array
553:
554: Level: advanced
556: .seealso VecPackDestroy(), VecPackGather(), VecPackAddDA(), VecPackCreateGlobalVector(),
557: VecPackScatter(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
558: VecPackGetLocalVectors(), VecPackRestoreLocalVectors()
560: @*/
561: PetscErrorCode VecPackAddArray(VecPack packer,PetscMPIInt orank,PetscInt n)
562: {
563: struct VecPackLink *mine,*next = packer->next;
564: PetscErrorCode ierr;
565: PetscMPIInt rank;
568: if (packer->globalvector) {
569: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot add an array once you have called VecPackCreateGlobalVector()");
570: }
571: #if defined(PETSC_USE_DEBUG)
572: {
573: PetscMPIInt orankmax;
574: MPI_Allreduce(&orank,&orankmax,1,MPI_INT,MPI_MAX,packer->comm);
575: if (orank != orankmax) SETERRQ2(PETSC_ERR_ARG_INCOMP,"orank %d must be equal on all processes, another process has value %d",orank,orankmax);
576: }
577: #endif
579: MPI_Comm_rank(packer->comm,&rank);
580: /* create new link */
581: PetscNew(struct VecPackLink,&mine);
582: mine->n = n;
583: mine->rank = orank;
584: mine->da = PETSC_NULL;
585: mine->type = VECPACK_ARRAY;
586: mine->next = PETSC_NULL;
587: if (rank == mine->rank) packer->n += n;
589: /* add to end of list */
590: if (!next) {
591: packer->next = mine;
592: } else {
593: while (next->next) next = next->next;
594: next->next = mine;
595: }
596: packer->nredundant++;
597: return(0);
598: }
602: /*@C
603: VecPackAddDA - adds a DA vector to a VecPack
605: Collective on VecPack
607: Input Parameter:
608: + packer - the packer object
609: - da - the DA object
610:
611: Level: advanced
613: .seealso VecPackDestroy(), VecPackGather(), VecPackAddDA(), VecPackCreateGlobalVector(),
614: VecPackScatter(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
615: VecPackGetLocalVectors(), VecPackRestoreLocalVectors()
617: @*/
618: PetscErrorCode VecPackAddDA(VecPack packer,DA da)
619: {
620: PetscErrorCode ierr;
621: PetscInt n;
622: struct VecPackLink *mine,*next = packer->next;
623: Vec global;
626: if (packer->globalvector) {
627: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DA once you have called VecPackCreateGlobalVector()");
628: }
630: /* create new link */
631: PetscNew(struct VecPackLink,&mine);
632: PetscObjectReference((PetscObject)da);
633: DAGetGlobalVector(da,&global);
634: VecGetLocalSize(global,&n);
635: DARestoreGlobalVector(da,&global);
636: mine->n = n;
637: mine->da = da;
638: mine->type = VECPACK_DA;
639: mine->next = PETSC_NULL;
640: packer->n += n;
642: /* add to end of list */
643: if (!next) {
644: packer->next = mine;
645: } else {
646: while (next->next) next = next->next;
647: next->next = mine;
648: }
649: packer->nDA++;
650: return(0);
651: }
655: /*@C
656: VecPackCreateGlobalVector - Creates a vector of the correct size to be gathered into
657: by the packer.
659: Collective on VecPack
661: Input Parameter:
662: . packer - the packer object
664: Output Parameters:
665: . gvec - the global vector
667: Level: advanced
669: Notes: Once this has been created you cannot add additional arrays or vectors to be packed.
671: This also serves as a VecPackFinalize(), perhaps that should be made a seperate
672: function?
674: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackScatter(),
675: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
676: VecPackGetLocalVectors(), VecPackRestoreLocalVectors()
678: @*/
679: PetscErrorCode VecPackCreateGlobalVector(VecPack packer,Vec *gvec)
680: {
681: PetscErrorCode ierr;
682: PetscInt nprev = 0;
683: PetscMPIInt rank,size;
684: struct VecPackLink *next = packer->next;
687: if (packer->globalvector) {
688: VecDuplicate(packer->globalvector,gvec);
689: } else {
690: VecCreateMPI(packer->comm,packer->n,PETSC_DETERMINE,gvec);
691: PetscObjectReference((PetscObject)*gvec);
692: packer->globalvector = *gvec;
694: VecGetSize(*gvec,&packer->N);
695: VecGetOwnershipRange(*gvec,&packer->rstart,PETSC_NULL);
696:
697: /* now set the rstart for each linked array/vector */
698: MPI_Comm_rank(packer->comm,&rank);
699: MPI_Comm_size(packer->comm,&size);
700: while (next) {
701: next->rstart = nprev;
702: if ((rank == next->rank) || next->type != VECPACK_ARRAY) nprev += next->n;
703: next->grstart = packer->rstart + next->rstart;
704: if (next->type == VECPACK_ARRAY) {
705: MPI_Bcast(&next->grstart,1,MPIU_INT,next->rank,packer->comm);
706: } else {
707: PetscMalloc(size*sizeof(PetscInt),&next->grstarts);
708: MPI_Allgather(&next->grstart,1,MPIU_INT,next->grstarts,1,MPIU_INT,packer->comm);
709: }
710: next = next->next;
711: }
712: VecSetOperation(*gvec,VECOP_VIEW,(void(*)(void))VecView_DADADA);
713: }
714: return(0);
715: }
719: /*@C
720: VecPackGetGlobalIndices - Gets the global indices for all the entries in the packed
721: vectors.
723: Collective on VecPack
725: Input Parameter:
726: . packer - the packer object
728: Output Parameters:
729: . idx - the individual indices for each packed vector/array. Note that this includes
730: all the ghost points that individual ghosted DA's may have.
731:
732: Level: advanced
734: Notes:
735: The idx parameters should be freed by the calling routine with PetscFree()
737: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
738: VecPackGather(), VecPackCreate(), VecPackGetAccess(),
739: VecPackGetLocalVectors(), VecPackRestoreLocalVectors()
741: @*/
742: PetscErrorCode VecPackGetGlobalIndices(VecPack packer,...)
743: {
744: va_list Argp;
745: PetscErrorCode ierr;
746: PetscInt i,**idx,n;
747: struct VecPackLink *next = packer->next;
748: Vec global,dglobal;
749: PF pf;
750: PetscScalar *array;
751: PetscMPIInt rank;
754: VecPackCreateGlobalVector(packer,&global);
755: MPI_Comm_rank(packer->comm,&rank);
757: /* put 0 to N-1 into the global vector */
758: PFCreate(PETSC_COMM_WORLD,1,1,&pf);
759: PFSetType(pf,PFIDENTITY,PETSC_NULL);
760: PFApplyVec(pf,PETSC_NULL,global);
761: PFDestroy(pf);
763: /* loop over packed objects, handling one at at time */
764: va_start(Argp,packer);
765: while (next) {
766: idx = va_arg(Argp, PetscInt**);
768: if (next->type == VECPACK_ARRAY) {
769:
770: PetscMalloc(next->n*sizeof(PetscInt),idx);
771: if (rank == next->rank) {
772: VecGetArray(global,&array);
773: array += next->rstart;
774: for (i=0; i<next->n; i++) (*idx)[i] = (PetscInt)PetscRealPart(array[i]);
775: array -= next->rstart;
776: VecRestoreArray(global,&array);
777: }
778: MPI_Bcast(*idx,next->n,MPIU_INT,next->rank,packer->comm);
780: } else if (next->type == VECPACK_DA) {
781: Vec local;
783: DACreateLocalVector(next->da,&local);
784: VecGetArray(global,&array);
785: array += next->rstart;
786: DAGetGlobalVector(next->da,&dglobal);
787: VecPlaceArray(dglobal,array);
788: DAGlobalToLocalBegin(next->da,dglobal,INSERT_VALUES,local);
789: DAGlobalToLocalEnd(next->da,dglobal,INSERT_VALUES,local);
790: array -= next->rstart;
791: VecRestoreArray(global,&array);
792: VecResetArray(dglobal);
793: DARestoreGlobalVector(next->da,&dglobal);
795: VecGetArray(local,&array);
796: VecGetSize(local,&n);
797: PetscMalloc(n*sizeof(PetscInt),idx);
798: for (i=0; i<n; i++) (*idx)[i] = (PetscInt)PetscRealPart(array[i]);
799: VecRestoreArray(local,&array);
800: VecDestroy(local);
802: } else {
803: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
804: }
805: next = next->next;
806: }
807: va_end(Argp);
808: VecDestroy(global);
809: return(0);
810: }
812: /* -------------------------------------------------------------------------------------*/
815: PetscErrorCode VecPackGetLocalVectors_Array(VecPack packer,struct VecPackLink *mine,PetscScalar **array)
816: {
819: PetscMalloc(mine->n*sizeof(PetscScalar),array);
820: return(0);
821: }
825: PetscErrorCode VecPackGetLocalVectors_DA(VecPack packer,struct VecPackLink *mine,Vec *local)
826: {
829: DAGetLocalVector(mine->da,local);
830: return(0);
831: }
835: PetscErrorCode VecPackRestoreLocalVectors_Array(VecPack packer,struct VecPackLink *mine,PetscScalar **array)
836: {
839: PetscFree(*array);
840: return(0);
841: }
845: PetscErrorCode VecPackRestoreLocalVectors_DA(VecPack packer,struct VecPackLink *mine,Vec *local)
846: {
849: DARestoreLocalVector(mine->da,local);
850: return(0);
851: }
855: /*@C
856: VecPackGetLocalVectors - Gets local vectors and arrays for each part of a VecPack.'
857: Use VecPakcRestoreLocalVectors() to return them.
859: Collective on VecPack
861: Input Parameter:
862: . packer - the packer object
863:
864: Output Parameter:
865: . ... - the individual sequential objects (arrays or vectors)
866:
867: Level: advanced
869: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
870: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
871: VecPackRestoreLocalVectors()
873: @*/
874: PetscErrorCode VecPackGetLocalVectors(VecPack packer,...)
875: {
876: va_list Argp;
877: PetscErrorCode ierr;
878: struct VecPackLink *next = packer->next;
882: /* loop over packed objects, handling one at at time */
883: va_start(Argp,packer);
884: while (next) {
885: if (next->type == VECPACK_ARRAY) {
886: PetscScalar **array;
887: array = va_arg(Argp, PetscScalar**);
888: VecPackGetLocalVectors_Array(packer,next,array);
889: } else if (next->type == VECPACK_DA) {
890: Vec *vec;
891: vec = va_arg(Argp, Vec*);
892: VecPackGetLocalVectors_DA(packer,next,vec);
893: } else {
894: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
895: }
896: next = next->next;
897: }
898: va_end(Argp);
899: return(0);
900: }
904: /*@C
905: VecPackRestoreLocalVectors - Restores local vectors and arrays for each part of a VecPack.'
906: Use VecPakcRestoreLocalVectors() to return them.
908: Collective on VecPack
910: Input Parameter:
911: . packer - the packer object
912:
913: Output Parameter:
914: . ... - the individual sequential objects (arrays or vectors)
915:
916: Level: advanced
918: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
919: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
920: VecPackGetLocalVectors()
922: @*/
923: PetscErrorCode VecPackRestoreLocalVectors(VecPack packer,...)
924: {
925: va_list Argp;
926: PetscErrorCode ierr;
927: struct VecPackLink *next = packer->next;
931: /* loop over packed objects, handling one at at time */
932: va_start(Argp,packer);
933: while (next) {
934: if (next->type == VECPACK_ARRAY) {
935: PetscScalar **array;
936: array = va_arg(Argp, PetscScalar**);
937: VecPackRestoreLocalVectors_Array(packer,next,array);
938: } else if (next->type == VECPACK_DA) {
939: Vec *vec;
940: vec = va_arg(Argp, Vec*);
941: VecPackRestoreLocalVectors_DA(packer,next,vec);
942: } else {
943: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
944: }
945: next = next->next;
946: }
947: va_end(Argp);
948: return(0);
949: }
951: /* -------------------------------------------------------------------------------------*/
954: PetscErrorCode VecPackGetEntries_Array(VecPack packer,struct VecPackLink *mine,PetscInt *n)
955: {
957: if (n) *n = mine->n;
958: return(0);
959: }
963: PetscErrorCode VecPackGetEntries_DA(VecPack packer,struct VecPackLink *mine,DA *da)
964: {
966: if (da) *da = mine->da;
967: return(0);
968: }
972: /*@C
973: VecPackGetEntries - Gets the DA, redundant size, etc for each entry in a VecPack.
974: Use VecPackRestoreEntries() to return them.
976: Collective on VecPack
978: Input Parameter:
979: . packer - the packer object
980:
981: Output Parameter:
982: . ... - the individual entries, DAs or integer sizes)
983:
984: Level: advanced
986: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
987: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
988: VecPackRestoreLocalVectors(), VecPackGetLocalVectors(), VecPackRestoreEntries(),
989: VecPackGetLocalVectors(), VecPackRestoreLocalVectors()
991: @*/
992: PetscErrorCode VecPackGetEntries(VecPack packer,...)
993: {
994: va_list Argp;
995: PetscErrorCode ierr;
996: struct VecPackLink *next = packer->next;
1000: /* loop over packed objects, handling one at at time */
1001: va_start(Argp,packer);
1002: while (next) {
1003: if (next->type == VECPACK_ARRAY) {
1004: PetscInt *n;
1005: n = va_arg(Argp, PetscInt*);
1006: VecPackGetEntries_Array(packer,next,n);
1007: } else if (next->type == VECPACK_DA) {
1008: DA *da;
1009: da = va_arg(Argp, DA*);
1010: VecPackGetEntries_DA(packer,next,da);
1011: } else {
1012: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1013: }
1014: next = next->next;
1015: }
1016: va_end(Argp);
1017: return(0);
1018: }
1022: /*@C
1023: VecPackRefine - Refines a VecPack by refining all of its DAs
1025: Collective on VecPack
1027: Input Parameters:
1028: + packer - the packer object
1029: - comm - communicator to contain the new DM object, usually PETSC_NULL
1031: Output Parameter:
1032: . fine - new packer
1033:
1034: Level: advanced
1036: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
1037: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
1038: VecPackGetLocalVectors(), VecPackRestoreLocalVectors()
1040: @*/
1041: PetscErrorCode VecPackRefine(VecPack packer,MPI_Comm comm,VecPack *fine)
1042: {
1043: PetscErrorCode ierr;
1044: struct VecPackLink *next = packer->next;
1045: DA da;
1048: VecPackCreate(comm,fine);
1050: /* loop over packed objects, handling one at at time */
1051: while (next) {
1052: if (next->type == VECPACK_ARRAY) {
1053: VecPackAddArray(*fine,next->rank,next->n);
1054: } else if (next->type == VECPACK_DA) {
1055: DARefine(next->da,comm,&da);
1056: VecPackAddDA(*fine,da);
1057: PetscObjectDereference((PetscObject)da);
1058: } else {
1059: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1060: }
1061: next = next->next;
1062: }
1063: return(0);
1064: }
1066: #include petscmat.h
1068: struct MatPackLink {
1069: Mat A;
1070: struct MatPackLink *next;
1071: };
1073: struct MatPack {
1074: VecPack right,left;
1075: struct MatPackLink *next;
1076: };
1080: PetscErrorCode MatMultBoth_Shell_Pack(Mat A,Vec x,Vec y,PetscTruth add)
1081: {
1082: struct MatPack *mpack;
1083: struct VecPackLink *xnext,*ynext;
1084: struct MatPackLink *anext;
1085: PetscScalar *xarray,*yarray;
1086: PetscErrorCode ierr;
1087: PetscInt i;
1088: Vec xglobal,yglobal;
1089: PetscMPIInt rank;
1092: MatShellGetContext(A,(void**)&mpack);
1093: MPI_Comm_rank(mpack->right->comm,&rank);
1094: xnext = mpack->right->next;
1095: ynext = mpack->left->next;
1096: anext = mpack->next;
1098: while (xnext) {
1099: if (xnext->type == VECPACK_ARRAY) {
1100: if (rank == xnext->rank) {
1101: VecGetArray(x,&xarray);
1102: VecGetArray(y,&yarray);
1103: if (add) {
1104: for (i=0; i<xnext->n; i++) {
1105: yarray[ynext->rstart+i] += xarray[xnext->rstart+i];
1106: }
1107: } else {
1108: PetscMemcpy(yarray+ynext->rstart,xarray+xnext->rstart,xnext->n*sizeof(PetscScalar));
1109: }
1110: VecRestoreArray(x,&xarray);
1111: VecRestoreArray(y,&yarray);
1112: }
1113: } else if (xnext->type == VECPACK_DA) {
1114: VecGetArray(x,&xarray);
1115: VecGetArray(y,&yarray);
1116: DAGetGlobalVector(xnext->da,&xglobal);
1117: DAGetGlobalVector(ynext->da,&yglobal);
1118: VecPlaceArray(xglobal,xarray+xnext->rstart);
1119: VecPlaceArray(yglobal,yarray+ynext->rstart);
1120: if (add) {
1121: MatMultAdd(anext->A,xglobal,yglobal,yglobal);
1122: } else {
1123: MatMult(anext->A,xglobal,yglobal);
1124: }
1125: VecRestoreArray(x,&xarray);
1126: VecRestoreArray(y,&yarray);
1127: VecResetArray(xglobal);
1128: VecResetArray(yglobal);
1129: DARestoreGlobalVector(xnext->da,&xglobal);
1130: DARestoreGlobalVector(ynext->da,&yglobal);
1131: anext = anext->next;
1132: } else {
1133: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1134: }
1135: xnext = xnext->next;
1136: ynext = ynext->next;
1137: }
1138: return(0);
1139: }
1143: PetscErrorCode MatMultAdd_Shell_Pack(Mat A,Vec x,Vec y,Vec z)
1144: {
1147: if (z != y) SETERRQ(PETSC_ERR_SUP,"Handles y == z only");
1148: MatMultBoth_Shell_Pack(A,x,y,PETSC_TRUE);
1149: return(0);
1150: }
1154: PetscErrorCode MatMult_Shell_Pack(Mat A,Vec x,Vec y)
1155: {
1158: MatMultBoth_Shell_Pack(A,x,y,PETSC_FALSE);
1159: return(0);
1160: }
1164: PetscErrorCode MatMultTranspose_Shell_Pack(Mat A,Vec x,Vec y)
1165: {
1166: struct MatPack *mpack;
1167: struct VecPackLink *xnext,*ynext;
1168: struct MatPackLink *anext;
1169: PetscScalar *xarray,*yarray;
1170: PetscErrorCode ierr;
1171: Vec xglobal,yglobal;
1172: PetscMPIInt rank;
1175: MatShellGetContext(A,(void**)&mpack);
1176: MPI_Comm_rank(mpack->right->comm,&rank);
1177: xnext = mpack->left->next;
1178: ynext = mpack->right->next;
1179: anext = mpack->next;
1181: while (xnext) {
1182: if (xnext->type == VECPACK_ARRAY) {
1183: if (rank == ynext->rank) {
1184: VecGetArray(x,&xarray);
1185: VecGetArray(y,&yarray);
1186: PetscMemcpy(yarray+ynext->rstart,xarray+xnext->rstart,xnext->n*sizeof(PetscScalar));
1187: VecRestoreArray(x,&xarray);
1188: VecRestoreArray(y,&yarray);
1189: }
1190: } else if (xnext->type == VECPACK_DA) {
1191: VecGetArray(x,&xarray);
1192: VecGetArray(y,&yarray);
1193: DAGetGlobalVector(xnext->da,&xglobal);
1194: DAGetGlobalVector(ynext->da,&yglobal);
1195: VecPlaceArray(xglobal,xarray+xnext->rstart);
1196: VecPlaceArray(yglobal,yarray+ynext->rstart);
1197: MatMultTranspose(anext->A,xglobal,yglobal);
1198: VecRestoreArray(x,&xarray);
1199: VecRestoreArray(y,&yarray);
1200: VecResetArray(xglobal);
1201: VecResetArray(yglobal);
1202: DARestoreGlobalVector(xnext->da,&xglobal);
1203: DARestoreGlobalVector(ynext->da,&yglobal);
1204: anext = anext->next;
1205: } else {
1206: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1207: }
1208: xnext = xnext->next;
1209: ynext = ynext->next;
1210: }
1211: return(0);
1212: }
1216: PetscErrorCode MatDestroy_Shell_Pack(Mat A)
1217: {
1218: struct MatPack *mpack;
1219: struct MatPackLink *anext,*oldanext;
1220: PetscErrorCode ierr;
1223: MatShellGetContext(A,(void**)&mpack);
1224: anext = mpack->next;
1226: while (anext) {
1227: MatDestroy(anext->A);
1228: oldanext = anext;
1229: anext = anext->next;
1230: PetscFree(oldanext);
1231: }
1232: PetscFree(mpack);
1233: PetscObjectChangeTypeName((PetscObject)A,0);
1234: return(0);
1235: }
1239: /*@C
1240: VecPackGetInterpolation - GetInterpolations a VecPack by refining all of its DAs
1242: Collective on VecPack
1244: Input Parameters:
1245: + coarse - coarse grid packer
1246: - fine - fine grid packer
1248: Output Parameter:
1249: + A - interpolation matrix
1250: - v - scaling vector
1251:
1252: Level: advanced
1254: .seealso VecPackDestroy(), VecPackAddArray(), VecPackAddDA(), VecPackCreateGlobalVector(),
1255: VecPackGather(), VecPackCreate(), VecPackGetGlobalIndices(), VecPackGetAccess(),
1256: VecPackGetLocalVectors(), VecPackRestoreLocalVectors()
1258: @*/
1259: PetscErrorCode VecPackGetInterpolation(VecPack coarse,VecPack fine,Mat *A,Vec *v)
1260: {
1261: PetscErrorCode ierr;
1262: PetscInt m,n,M,N;
1263: struct VecPackLink *nextc = coarse->next;
1264: struct VecPackLink *nextf = fine->next;
1265: struct MatPackLink *nextmat,*pnextmat = 0;
1266: struct MatPack *mpack;
1267: Vec gcoarse,gfine;
1270: /* use global vectors only for determining matrix layout */
1271: VecPackCreateGlobalVector(coarse,&gcoarse);
1272: VecPackCreateGlobalVector(fine,&gfine);
1273: VecGetLocalSize(gcoarse,&n);
1274: VecGetLocalSize(gfine,&m);
1275: VecGetSize(gcoarse,&N);
1276: VecGetSize(gfine,&M);
1277: VecDestroy(gcoarse);
1278: VecDestroy(gfine);
1280: PetscNew(struct MatPack,&mpack);
1281: mpack->right = coarse;
1282: mpack->left = fine;
1283: MatCreate(fine->comm,A);
1284: MatSetSizes(*A,m,n,M,N);
1285: MatSetType(*A,MATSHELL);
1286: MatShellSetContext(*A,mpack);
1287: MatShellSetOperation(*A,MATOP_MULT,(void(*)(void))MatMult_Shell_Pack);
1288: MatShellSetOperation(*A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Shell_Pack);
1289: MatShellSetOperation(*A,MATOP_MULT_ADD,(void(*)(void))MatMultAdd_Shell_Pack);
1290: MatShellSetOperation(*A,MATOP_DESTROY,(void(*)(void))MatDestroy_Shell_Pack);
1292: /* loop over packed objects, handling one at at time */
1293: while (nextc) {
1294: if (nextc->type != nextf->type) SETERRQ(PETSC_ERR_ARG_INCOMP,"Two VecPack have different layout");
1296: if (nextc->type == VECPACK_ARRAY) {
1297: ;
1298: } else if (nextc->type == VECPACK_DA) {
1299: PetscNew(struct MatPackLink,&nextmat);
1300: nextmat->next = 0;
1301: if (pnextmat) {
1302: pnextmat->next = nextmat;
1303: pnextmat = nextmat;
1304: } else {
1305: pnextmat = nextmat;
1306: mpack->next = nextmat;
1307: }
1308: DAGetInterpolation(nextc->da,nextf->da,&nextmat->A,PETSC_NULL);
1309: } else {
1310: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1311: }
1312: nextc = nextc->next;
1313: nextf = nextf->next;
1314: }
1315: return(0);
1316: }
1320: /*@C
1321: VecPackGetMatrix - Creates a matrix with the correct parallel layout and nonzero structure required for
1322: computing the Jacobian on a function defined using the stencils set in the DA's and coupling in the array variables
1324: Collective on DA
1326: Input Parameter:
1327: + da - the distributed array
1328: - mtype - Supported types are MATSEQAIJ, MATMPIAIJ
1330: Output Parameters:
1331: . J - matrix with the correct nonzero structure
1332: (obviously without the correct Jacobian values)
1334: Level: advanced
1336: Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
1337: do not need to do it yourself.
1340: .seealso DAGetMatrix(), VecPackCreate()
1342: @*/
1343: PetscErrorCode VecPackGetMatrix(VecPack packer, MatType mtype,Mat *J)
1344: {
1345: PetscErrorCode ierr;
1346: struct VecPackLink *next = packer->next;
1347: Vec gvec;
1348: PetscInt m,*dnz,*onz,i,j,mA;
1349: Mat Atmp;
1350: PetscMPIInt rank;
1351: PetscScalar zero = 0.0;
1354: /* use global vector to determine layout needed for matrix */
1355: VecPackCreateGlobalVector(packer,&gvec);
1356: VecGetLocalSize(gvec,&m);
1357: VecDestroy(gvec);
1359: MPI_Comm_rank(packer->comm,&rank);
1360: MatCreate(packer->comm,J);
1361: MatSetSizes(*J,m,m,PETSC_DETERMINE,PETSC_DETERMINE);
1362: MatSetType(*J,mtype);
1364: MatPreallocateInitialize(packer->comm,m,m,dnz,onz);
1365: /* loop over packed objects, handling one at at time */
1366: while (next) {
1367: if (next->type == VECPACK_ARRAY) {
1368: if (rank == next->rank) { /* zero the entire row */
1369: for (j=packer->rstart+next->rstart; j<packer->rstart+next->rstart+next->n; j++) {
1370: for (i=0; i<packer->N; i++) {
1371: MatPreallocateSet(j,1,&i,dnz,onz);
1372: }
1373: }
1374: }
1375: for (j=next->grstart; j<next->grstart+next->n; j++) {
1376: for (i=packer->rstart; i<packer->rstart+m; i++) { /* zero the entire local column */
1377: if (j != i) { /* don't count diagonal twice */
1378: MatPreallocateSet(i,1,&j,dnz,onz);
1379: }
1380: }
1381: }
1382: } else if (next->type == VECPACK_DA) {
1383: PetscInt nc,rstart,*ccols,maxnc;
1384: const PetscInt *cols,*rstarts;
1385: PetscMPIInt proc;
1387: DAGetMatrix(next->da,mtype,&Atmp);
1388: MatGetOwnershipRange(Atmp,&rstart,PETSC_NULL);
1389: MatGetOwnershipRanges(Atmp,&rstarts);
1390: MatGetLocalSize(Atmp,&mA,PETSC_NULL);
1392: maxnc = 0;
1393: for (i=0; i<mA; i++) {
1394: MatGetRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);
1395: MatRestoreRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);
1396: maxnc = PetscMax(nc,maxnc);
1397: }
1398: PetscMalloc(maxnc*sizeof(PetscInt),&ccols);
1399: for (i=0; i<mA; i++) {
1400: MatGetRow(Atmp,rstart+i,&nc,&cols,PETSC_NULL);
1401: /* remap the columns taking into how much they are shifted on each process */
1402: for (j=0; j<nc; j++) {
1403: proc = 0;
1404: while (cols[j] >= rstarts[proc+1]) proc++;
1405: ccols[j] = cols[j] + next->grstarts[proc] - rstarts[proc];
1406: }
1407: MatPreallocateSet(packer->rstart+next->rstart+i,nc,ccols,dnz,onz);
1408: MatRestoreRow(Atmp,rstart+i,&nc,&cols,PETSC_NULL);
1409: }
1410: PetscFree(ccols);
1411: } else {
1412: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1413: }
1414: next = next->next;
1415: }
1416: MatMPIAIJSetPreallocation(*J,0,dnz,0,onz);
1417: MatSeqAIJSetPreallocation(*J,0,dnz);
1418: MatPreallocateFinalize(dnz,onz);
1420: next = packer->next;
1421: while (next) {
1422: if (next->type == VECPACK_ARRAY) {
1423: if (rank == next->rank) {
1424: for (j=packer->rstart+next->rstart; j<packer->rstart+next->rstart+next->n; j++) {
1425: for (i=0; i<packer->N; i++) {
1426: MatSetValues(*J,1,&j,1,&i,&zero,INSERT_VALUES);
1427: }
1428: }
1429: }
1430: for (j=next->grstart; j<next->grstart+next->n; j++) {
1431: for (i=packer->rstart; i<packer->rstart+m; i++) {
1432: MatSetValues(*J,1,&i,1,&j,&zero,INSERT_VALUES);
1433: }
1434: }
1435: } else if (next->type == VECPACK_DA) {
1436: PetscInt nc,rstart,row,maxnc,*ccols;
1437: const PetscInt *cols,*rstarts;
1438: const PetscScalar *values;
1439: PetscMPIInt proc;
1441: DAGetMatrix(next->da,mtype,&Atmp);
1442: MatGetOwnershipRange(Atmp,&rstart,PETSC_NULL);
1443: MatGetOwnershipRanges(Atmp,&rstarts);
1444: MatGetLocalSize(Atmp,&mA,PETSC_NULL);
1445: maxnc = 0;
1446: for (i=0; i<mA; i++) {
1447: MatGetRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);
1448: MatRestoreRow(Atmp,rstart+i,&nc,PETSC_NULL,PETSC_NULL);
1449: maxnc = PetscMax(nc,maxnc);
1450: }
1451: PetscMalloc(maxnc*sizeof(PetscInt),&ccols);
1452: for (i=0; i<mA; i++) {
1453: MatGetRow(Atmp,rstart+i,&nc,(const PetscInt **)&cols,&values);
1454: for (j=0; j<nc; j++) {
1455: proc = 0;
1456: while (cols[j] >= rstarts[proc+1]) proc++;
1457: ccols[j] = cols[j] + next->grstarts[proc] - rstarts[proc];
1458: }
1459: row = packer->rstart+next->rstart+i;
1460: MatSetValues(*J,1,&row,nc,ccols,values,INSERT_VALUES);
1461: MatRestoreRow(Atmp,rstart+i,&nc,(const PetscInt **)&cols,&values);
1462: }
1463: PetscFree(ccols);
1464: } else {
1465: SETERRQ(PETSC_ERR_SUP,"Cannot handle that object type yet");
1466: }
1467: next = next->next;
1468: }
1469: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
1470: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
1471: return(0);
1472: }
1476: /*@
1477: VecPackGetColoring - Gets the coloring required for computing the Jacobian via
1478: finite differences on a function defined using a VecPack "grid"
1480: Collective on DA
1482: Input Parameter:
1483: + vecpack - the VecPack object
1484: - ctype - IS_COLORING_GLOBAL or IS_COLORING_GHOSTED
1486: Output Parameters:
1487: . coloring - matrix coloring for use in computing Jacobians (or PETSC_NULL if not needed)
1489: Level: advanced
1491: Notes: This currentlu uses one color per column so is very slow.
1493: Notes: These compute the graph coloring of the graph of A^{T}A. The coloring used
1494: for efficient (parallel or thread based) triangular solves etc is NOT yet
1495: available.
1498: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), ISColoringType, ISColoring, DAGetColoring()
1500: @*/
1501: PetscErrorCode VecPackGetColoring(VecPack vecpack,ISColoringType ctype,ISColoring *coloring)
1502: {
1503: PetscErrorCode ierr;
1504: PetscInt n,i;
1505: ISColoringValue *colors;
1508: if (ctype == IS_COLORING_GHOSTED) {
1509: SETERRQ(PETSC_ERR_SUP,"Lazy Barry");
1510: } else if (ctype == IS_COLORING_GLOBAL) {
1511: n = vecpack->n;
1512: } else SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1514: PetscMalloc(n*sizeof(ISColoringValue),&colors); /* freed in ISColoringDestroy() */
1515: for (i=0; i<n; i++) {
1516: colors[i] = (ISColoringValue)(vecpack->rstart + i);
1517: }
1518: ISColoringCreate(vecpack->comm,vecpack->N,n,colors,coloring);
1519: return(0);
1520: }