Actual source code: isutil.c
1: #include "tao.h" /*I "tao.h" I*/
2: #include "private/matimpl.h"
3: #include "src/matrix/submatfree.h"
4: #include "tao_util.h" /*I "tao_util.h" I*/
9: /*@
10: VecWhichEqual - Creates an index set containing the indices
11: where the vectors Vec1 and Vec2 have identical elements.
12:
13: Collective on S
15: Input Parameters:
16: . Vec1, Vec2 - the two vectors to compare
18: OutputParameter:
19: . S - The index set containing the indices i where vec1[i] == vec2[i]
21: Level: advanced
22: @*/
23: PetscErrorCode VecWhichEqual(Vec Vec1, Vec Vec2, IS * S)
24: {
25: PetscErrorCode ierr;
26: PetscInt i,n_same=0;
27: PetscInt n,low,high,low2,high2;
28: PetscInt *same;
29: PetscReal *v1,*v2;
30: MPI_Comm comm;
38: VecGetOwnershipRange(Vec1, &low, &high);
39: VecGetOwnershipRange(Vec2, &low2, &high2);
41: if ( low != low2 || high != high2 )
42: SETERRQ(PETSC_COMM_SELF,1,"Vectors must be identically loaded over processors");
44: VecGetLocalSize(Vec1,&n);
46: if (n>0){
47:
48: if (Vec1 == Vec2){
49: VecGetArray(Vec1,&v1);
50: v2=v1;
51: } else {
52: VecGetArray(Vec1,&v1);
53: VecGetArray(Vec2,&v2);
54: }
56: PetscMalloc( n*sizeof(PetscInt),&same );
57:
58: for (i=0; i<n; i++){
59: if (v1[i] == v2[i]) {same[n_same]=low+i; n_same++;}
60: }
61:
62: if (Vec1 == Vec2){
63: VecRestoreArray(Vec1,&v1);
64: } else {
65: VecRestoreArray(Vec1,&v1);
66: VecRestoreArray(Vec2,&v2);
67: }
69: } else {
71: n_same = 0; same=NULL;
73: }
75: PetscObjectGetComm((PetscObject)Vec1,&comm);
76: ISCreateGeneral(comm,n_same,same,PETSC_COPY_VALUES,S);
78: if (same) {
79: PetscFree(same);
80: }
82: return(0);
83: }
87: /*@
88: VecWhichLessThan - Creates an index set containing the indices
89: where the vectors Vec1 < Vec2
90:
91: Collective on S
93: Input Parameters:
94: . Vec1, Vec2 - the two vectors to compare
96: OutputParameter:
97: . S - The index set containing the indices i where vec1[i] < vec2[i]
99: Level: advanced
100: @*/
101: PetscErrorCode VecWhichLessThan(Vec Vec1, Vec Vec2, IS * S)
102: {
103: int ierr;
104: PetscInt i;
105: PetscInt n,low,high,low2,high2,n_lt=0;
106: PetscInt *lt;
107: PetscReal *v1,*v2;
108: MPI_Comm comm;
115: VecGetOwnershipRange(Vec1, &low, &high);
116: VecGetOwnershipRange(Vec2, &low2, &high2);
118: if ( low != low2 || high != high2 )
119: SETERRQ(PETSC_COMM_SELF,1,"Vectors must be identically loaded over processors");
121: VecGetLocalSize(Vec1,&n);
123: if (n>0){
125: if (Vec1 == Vec2){
126: VecGetArray(Vec1,&v1);
127: v2=v1;
128: } else {
129: VecGetArray(Vec1,&v1);
130: VecGetArray(Vec2,&v2);
131: }
132: PetscMalloc( n*sizeof(PetscInt),< );
133:
134: for (i=0; i<n; i++){
135: if (v1[i] < v2[i]) {lt[n_lt]=high+i; n_lt++;}
136: }
138: if (Vec1 == Vec2){
139: VecRestoreArray(Vec1,&v1);
140: } else {
141: VecRestoreArray(Vec1,&v1);
142: VecRestoreArray(Vec2,&v2);
143: }
144:
145: } else {
146: n_lt=0; lt=NULL;
147: }
149: PetscObjectGetComm((PetscObject)Vec1,&comm);
150: ISCreateGeneral(comm,n_lt,lt,PETSC_COPY_VALUES,S);
152: if (lt) {
153: PetscFree(lt);
154: }
156: return(0);
157: }
161: /*@
162: VecWhichGreaterThan - Creates an index set containing the indices
163: where the vectors Vec1 > Vec2
164:
165: Collective on S
167: Input Parameters:
168: . Vec1, Vec2 - the two vectors to compare
170: OutputParameter:
171: . S - The index set containing the indices i where vec1[i] > vec2[i]
173: Level: advanced
174: @*/
175: PetscErrorCode VecWhichGreaterThan(Vec Vec1, Vec Vec2, IS * S)
176: {
177: int ierr;
178: PetscInt n,low,high,low2,high2,n_gt=0,i;
179: PetscInt *gt=NULL;
180: PetscReal *v1,*v2;
181: MPI_Comm comm;
188: VecGetOwnershipRange(Vec1, &low, &high);
189: VecGetOwnershipRange(Vec2, &low2, &high2);
191: if ( low != low2 || high != high2 )
192: SETERRQ(PETSC_COMM_SELF,1,"Vectors must be identically loaded over processors");
194: VecGetLocalSize(Vec1,&n);
196: if (n>0){
198: if (Vec1 == Vec2){
199: VecGetArray(Vec1,&v1);
200: v2=v1;
201: } else {
202: VecGetArray(Vec1,&v1);
203: VecGetArray(Vec2,&v2);
204: }
206: PetscMalloc( n*sizeof(PetscInt), > );
207:
208: for (i=0; i<n; i++){
209: if (v1[i] > v2[i]) {gt[n_gt]=low+i; n_gt++;}
210: }
212: if (Vec1 == Vec2){
213: VecRestoreArray(Vec1,&v1);
214: } else {
215: VecRestoreArray(Vec1,&v1);
216: VecRestoreArray(Vec2,&v2);
217: }
218:
219: } else{
220:
221: n_gt=0; gt=NULL;
223: }
225: PetscObjectGetComm((PetscObject)Vec1,&comm);
226: ISCreateGeneral(comm,n_gt,gt,PETSC_COPY_VALUES,S);
228: if (gt) {
229: PetscFree(gt);
230: }
231: return(0);
232: }
236: /*@
237: VecWhichBetween - Creates an index set containing the indices
238: where VecLow < V < VecHigh
239:
240: Collective on S
242: Input Parameters:
243: + VecLow - lower bound
244: . V - Vector to compare
245: - VecHigh - higher bound
247: OutputParameter:
248: . S - The index set containing the indices i where veclow[i] < v[i] < vechigh[i]
250: Level: advanced
251: @*/
252: PetscErrorCode VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS *S)
253: {
256: PetscInt n,low,high,low2,high2,low3,high3,n_vm=0;
257: PetscInt *vm,i;
258: PetscReal *v1,*v2,*vmiddle;
259: MPI_Comm comm;
264: VecGetOwnershipRange(VecLow, &low, &high);
265: VecGetOwnershipRange(VecHigh, &low2, &high2);
266: VecGetOwnershipRange(V, &low3, &high3);
268: if ( low!=low2 || high!=high2 || low!=low3 || high!=high3 )
269: SETERRQ(PETSC_COMM_SELF,1,"Vectors must be identically loaded over processors");
271: VecGetLocalSize(VecLow,&n);
273: if (n>0){
275: VecGetArray(VecLow,&v1);
276: if (VecLow != VecHigh){
277: VecGetArray(VecHigh,&v2);
278: } else {
279: v2=v1;
280: }
281: if ( V != VecLow && V != VecHigh){
282: VecGetArray(V,&vmiddle);
283: } else if ( V==VecLow ){
284: vmiddle=v1;
285: } else {
286: vmiddle =v2;
287: }
289: PetscMalloc( n*sizeof(PetscInt), &vm );
290:
291: for (i=0; i<n; i++){
292: if (v1[i] < vmiddle[i] && vmiddle[i] < v2[i]) {vm[n_vm]=low+i; n_vm++;}
293: }
295: VecRestoreArray(VecLow,&v1);
296: if (VecLow != VecHigh){
297: VecRestoreArray(VecHigh,&v2);
298: }
299: if ( V != VecLow && V != VecHigh){
300: VecRestoreArray(V,&vmiddle);
301: }
303: } else {
305: n_vm=0; vm=NULL;
307: }
309: PetscObjectGetComm((PetscObject)V,&comm);
310: ISCreateGeneral(comm,n_vm,vm,PETSC_COPY_VALUES,S);
312: if (vm) {
313: PetscFree(vm);
314: }
316: return(0);
317: }
322: /*@
323: VecWhichBetweenOrEqual - Creates an index set containing the indices
324: where VecLow <= V <= VecHigh
325:
326: Collective on S
328: Input Parameters:
329: + VecLow - lower bound
330: . V - Vector to compare
331: - VecHigh - higher bound
333: OutputParameter:
334: . S - The index set containing the indices i where veclow[i] <= v[i] <= vechigh[i]
336: Level: advanced
337: @*/
339: PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS * S)
340: {
342: PetscInt n,low,high,low2,high2,low3,high3,n_vm=0,i;
343: PetscInt *vm;
344: PetscReal *v1,*v2,*vmiddle;
345: MPI_Comm comm;
350: VecGetOwnershipRange(VecLow, &low, &high);
351: VecGetOwnershipRange(VecHigh, &low2, &high2);
352: VecGetOwnershipRange(V, &low3, &high3);
354: if ( low!=low2 || high!=high2 || low!=low3 || high!=high3 )
355: SETERRQ(PETSC_COMM_SELF,1,"Vectors must be identically loaded over processors");
357: VecGetLocalSize(VecLow,&n);
359: if (n>0){
361: VecGetArray(VecLow,&v1);
362: if (VecLow != VecHigh){
363: VecGetArray(VecHigh,&v2);
364: } else {
365: v2=v1;
366: }
367: if ( V != VecLow && V != VecHigh){
368: VecGetArray(V,&vmiddle);
369: } else if ( V==VecLow ){
370: vmiddle=v1;
371: } else {
372: vmiddle =v2;
373: }
375: PetscMalloc( n*sizeof(PetscInt), &vm );
376:
377: for (i=0; i<n; i++){
378: if (v1[i] <= vmiddle[i] && vmiddle[i] <= v2[i]) {vm[n_vm]=low+i; n_vm++;}
379: }
381: VecRestoreArray(VecLow,&v1);
382: if (VecLow != VecHigh){
383: VecRestoreArray(VecHigh,&v2);
384: }
385: if ( V != VecLow && V != VecHigh){
386: VecRestoreArray(V,&vmiddle);
387: }
389: } else {
391: n_vm=0; vm=NULL;
393: }
395: PetscObjectGetComm((PetscObject)V,&comm);
396: ISCreateGeneral(comm,n_vm,vm,PETSC_COPY_VALUES,S);
398: if (vm) {
399: PetscFree(vm);
400: }
402: return(0);
403: }
408: /*@
409: VecGetSubVec - Gets a subvector using the IS
411: Input Parameters:
412: + vfull - the full matrix
413: . is - the index set for the subvector
414: . reduced_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK, TAO_SUBSET_MATRIXFREE)
415: - maskvalue - the value to set the unused vector elements to (for TAO_SUBSET_MASK or TAO_SUBSET_MATRIXFREE)
418: Output Parameters:
419: . vreduced - the subvector
421: Note:
422: maskvalue should usually be 0.0, unless a pointwise divide will be used.
423: @*/
424: PetscErrorCode VecGetSubVec(Vec vfull, IS is, PetscInt reduced_type, PetscReal maskvalue, Vec *vreduced)
425: {
427: PetscInt nfull,nreduced,nreduced_local,rlow,rhigh,flow,fhigh;
428: PetscInt i,nlocal;
429: PetscReal *fv,*rv;
430: const PetscInt *s;
431: IS ident;
432: const VecType vtype;
433: VecScatter scatter;
434: MPI_Comm comm;
435:
439:
440: VecGetSize(vfull, &nfull);
441: ISGetSize(is, &nreduced);
443: if (nreduced == nfull) {
445: if (*vreduced == PETSC_NULL) {
446: VecDuplicate(vfull,vreduced);
447: }
448: VecDuplicate(vfull,vreduced);
449: VecCopy(vfull,*vreduced);
451: } else {
452:
453: switch (reduced_type) {
454: case TAO_SUBSET_SUBVEC:
455: VecGetType(vfull,&vtype);
456: VecGetOwnershipRange(vfull,&flow,&fhigh);
457: ISGetLocalSize(is,&nreduced_local);
458: PetscObjectGetComm((PetscObject)vfull,&comm);
459: if (*vreduced) {
460: VecDestroy(vreduced);
461: }
462: VecCreate(comm,vreduced);
463: VecSetType(*vreduced,vtype);
465:
466: VecSetSizes(*vreduced,nreduced_local,nreduced);
467:
468: VecGetOwnershipRange(*vreduced,&rlow,&rhigh);
469:
470: ISCreateStride(comm,nreduced_local,rlow,1,&ident);
471: VecScatterCreate(vfull,is,*vreduced,ident,&scatter);
472: VecScatterBegin(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);
473: VecScatterEnd(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);
474: VecScatterDestroy(&scatter);
475: ISDestroy(&ident);
476: break;
478: case TAO_SUBSET_MASK:
479: case TAO_SUBSET_MATRIXFREE:
480: /* vr[i] = vf[i] if i in is
481: vr[i] = 0 otherwise */
482: if (*vreduced == PETSC_NULL) {
483: VecDuplicate(vfull,vreduced);
484: }
486: VecSet(*vreduced,maskvalue);
487: ISGetLocalSize(is,&nlocal);
488: VecGetOwnershipRange(vfull,&flow,&fhigh);
489: VecGetArray(vfull,&fv);
490: VecGetArray(*vreduced,&rv);
491: ISGetIndices(is,&s);
492: if (nlocal > (fhigh-flow)) {
493: SETERRQ2(PETSC_COMM_WORLD,1,"IS local size %d > Vec local size %d",nlocal,fhigh-flow);
494: }
495: for (i=0;i<nlocal;i++) {
496: rv[s[i]-flow] = fv[s[i]-flow];
497: }
498: ISRestoreIndices(is,&s);
499: VecRestoreArray(vfull,&fv);
500: VecRestoreArray(*vreduced,&rv);
501: break;
502: }
503: }
504: return(0);
505:
506: }
510: /*@
511: VecReducedXPY - Adds a reduced vector to the appropriate elements of a full-space vector.
513: Input Parameters:
514: + vfull - the full-space vector
515: . vreduced - the reduced-space vector
516: - is - the index set for the reduced space
518: Output Parameters:
519: . vfull - the sum of the full-space vector and reduced-space vector
520: @*/
521: PetscErrorCode VecReducedXPY(Vec vfull, Vec vreduced, IS is)
522: {
523: VecScatter scatter;
524: IS ident;
525: PetscInt nfull,nreduced,rlow,rhigh;
526: MPI_Comm comm;
533: VecGetSize(vfull,&nfull);
534: VecGetSize(vreduced,&nreduced);
535:
536: if (nfull == nreduced) { /* Also takes care of masked vectors */
537: VecAXPY(vfull,1.0,vreduced);
538: } else {
539: PetscObjectGetComm((PetscObject)vfull,&comm);
540: VecGetOwnershipRange(vreduced,&rlow,&rhigh);
541: ISCreateStride(comm,rhigh-rlow,rlow,1,&ident);
542: VecScatterCreate(vreduced,ident,vfull,is,&scatter);
543: VecScatterBegin(scatter,vreduced,vfull,ADD_VALUES,SCATTER_FORWARD);
544: VecScatterEnd(scatter,vreduced,vfull,ADD_VALUES,SCATTER_FORWARD);
545: VecScatterDestroy(&scatter);
546: ISDestroy(&ident);
547: }
548:
549: return(0);
550: }
555: /*@
556: ISCreateComplement - Creates the complement of the the index set
558: Collective on IS
560: Input Parameter:
561: + S - a PETSc IS
562: - V - the reference vector space
564: Output Parameter:
565: . T - the complement of S
568: .seealso ISCreateGeneral()
570: Level: advanced
571: @*/
572: PetscErrorCode ISCreateComplement(IS S, Vec V, IS *T){
574: PetscInt i,nis,nloc,high,low,n=0;
575: const PetscInt *s;
576: PetscInt *tt,*ss;
577: MPI_Comm comm;
583: VecGetOwnershipRange(V,&low,&high);
584: VecGetLocalSize(V,&nloc);
585: ISGetLocalSize(S,&nis);
586: ISGetIndices(S, &s);
587: PetscMalloc( nloc*sizeof(PetscInt),&tt );
588: PetscMalloc( nloc*sizeof(PetscInt),&ss );
590: for (i=low; i<high; i++){ tt[i-low]=i; }
592: for (i=0; i<nis; i++){ tt[s[i]-low] = -2; }
593:
594: for (i=0; i<nloc; i++){
595: if (tt[i]>-1){ ss[n]=tt[i]; n++; }
596: }
598: ISRestoreIndices(S, &s);
599:
600: PetscObjectGetComm((PetscObject)S,&comm);
601: ISCreateGeneral(comm,n,ss,PETSC_COPY_VALUES,T);
602:
603: if (tt) {
604: PetscFree(tt);
605: }
606: if (ss) {
607: PetscFree(ss);
608: }
610: return(0);
611: }
615: /*@
616: VecISSetToConstant - Sets the elements of a vector, specified by an index set, to a constant
618: Input Parameter:
619: + S - a PETSc IS
620: . c - the constant
621: - V - a Vec
623: .seealso VecSet()
625: Level: advanced
626: @*/
627: PetscErrorCode VecISSetToConstant(IS S, PetscReal c, Vec V){
629: PetscInt nloc,low,high,i;
630: const PetscInt *s;
631: PetscReal *v;
638: VecGetOwnershipRange(V, &low, &high);
639: ISGetLocalSize(S,&nloc);
641: ISGetIndices(S, &s);
642: VecGetArray(V,&v);
643: for (i=0; i<nloc; i++){
644: v[s[i]-low] = c;
645: }
646:
647: ISRestoreIndices(S, &s);
648: VecRestoreArray(V,&v);
650: return(0);
652: }
656: /*@
657: MatGetSubMat - Gets a submatrix using the IS
659: Input Parameters:
660: + M - the full matrix (n x n)
661: . is - the index set for the submatrix (both row and column index sets need to be the same)
662: . v1 - work vector of dimension n, needed for TAO_SUBSET_MASK option
663: - subset_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK,
664: TAO_SUBSET_MATRIXFREE)
666: Output Parameters:
667: . Msub - the submatrix
668: @*/
669: PetscErrorCode MatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
670: {
672: IS iscomp;
673: PetscBool flg;
677: if (*Msub) {
678: MatDestroy(Msub);
679: }
680: switch (subset_type) {
681: case TAO_SUBSET_SUBVEC:
682: MatGetSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);
683: break;
685: case TAO_SUBSET_MASK:
686: /* Get Reduced Hessian
687: Msub[i,j] = M[i,j] if i,j in Free_Local or i==j
688: Msub[i,j] = 0 if i!=j and i or j not in Free_Local
689: */
690: PetscOptionsBool("-different_submatrix","use separate hessian matrix when computing submatrices","TaoSubsetType",PETSC_FALSE,&flg,PETSC_NULL);
691: if (flg == PETSC_TRUE) {
692: MatDuplicate(M, MAT_COPY_VALUES, Msub);
693: } else {
694: /* Act on hessian directly (default) */
695: PetscObjectReference((PetscObject)M);
696: *Msub = M;
697: }
698: /* Save the diagonal to temporary vector */
699: MatGetDiagonal(*Msub,v1);
700:
701: /* Zero out rows and columns */
702: ISCreateComplement(is,v1,&iscomp);
704: /* Use v1 instead of 0 here because of PETSc bug */
705: MatZeroRowsColumnsIS(*Msub,iscomp,1.0,v1,v1);
707: ISDestroy(&iscomp);
708: break;
709: case TAO_SUBSET_MATRIXFREE:
710: ISCreateComplement(is,v1,&iscomp);
711: MatCreateSubMatrixFree(M,iscomp,iscomp,Msub);
712: ISDestroy(&iscomp);
713: break;
714: }
715: return(0);
716: }