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),&lt ); 
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), &gt ); 
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: }