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:       VecDestroy(vreduced); 
446:       VecDuplicate(vfull,vreduced); 
447:       VecCopy(vfull,*vreduced); 

449:     } else { 
450:      
451:         switch (reduced_type) {
452:         case TAO_SUBSET_SUBVEC:
453:           VecGetType(vfull,&vtype); 
454:           VecGetOwnershipRange(vfull,&flow,&fhigh); 
455:           ISGetLocalSize(is,&nreduced_local); 
456:           PetscObjectGetComm((PetscObject)vfull,&comm); 
457:           if (*vreduced) {
458:             VecDestroy(vreduced); 
459:           }
460:           VecCreate(comm,vreduced); 
461:           VecSetType(*vreduced,vtype); 

463:           
464:           VecSetSizes(*vreduced,nreduced_local,nreduced); 
465:         
466:           VecGetOwnershipRange(*vreduced,&rlow,&rhigh); 
467:         
468:           ISCreateStride(comm,nreduced_local,rlow,1,&ident); 
469:           VecScatterCreate(vfull,is,*vreduced,ident,&scatter); 
470:           VecScatterBegin(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD); 
471:           VecScatterEnd(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD); 
472:           VecScatterDestroy(&scatter); 
473:           ISDestroy(&ident); 
474:           break;

476:         case TAO_SUBSET_MASK:
477:         case TAO_SUBSET_MATRIXFREE:
478:           /* vr[i] = vf[i]   if i in is
479:              vr[i] = 0       otherwise */
480:           if (*vreduced == PETSC_NULL) {
481:             VecDuplicate(vfull,vreduced); 
482:           }

484:           VecSet(*vreduced,maskvalue); 
485:           ISGetLocalSize(is,&nlocal); 
486:           VecGetOwnershipRange(vfull,&flow,&fhigh); 
487:           VecGetArray(vfull,&fv); 
488:           VecGetArray(*vreduced,&rv); 
489:           ISGetIndices(is,&s); 
490:           if (nlocal > (fhigh-flow)) {
491:             SETERRQ2(PETSC_COMM_WORLD,1,"IS local size %d > Vec local size %d",nlocal,fhigh-flow);
492:           }
493:           for (i=0;i<nlocal;i++) {
494:             rv[s[i]-flow] = fv[s[i]-flow];
495:           }
496:           ISRestoreIndices(is,&s); 
497:           VecRestoreArray(vfull,&fv); 
498:           VecRestoreArray(*vreduced,&rv); 
499:           break;
500:         }
501:     }
502:     return(0);
503:     
504: }

508: /*@ 
509:   VecReducedXPY - Adds a reduced vector to the appropriate elements of a full-space vector.

511:   Input Parameters:
512: + vfull - the full-space vector
513: . vreduced - the reduced-space vector
514: - is - the index set for the reduced space

516:   Output Parameters:
517: . vfull - the sum of the full-space vector and reduced-space vector
518: @*/
519: PetscErrorCode VecReducedXPY(Vec vfull, Vec vreduced, IS is)
520: {
521:     VecScatter scatter;
522:     IS ident;
523:     PetscInt nfull,nreduced,rlow,rhigh;
524:     MPI_Comm comm;

531:     VecGetSize(vfull,&nfull); 
532:     VecGetSize(vreduced,&nreduced); 
533:     
534:     if (nfull == nreduced) { /* Also takes care of masked vectors */
535:       VecAXPY(vfull,1.0,vreduced); 
536:     } else {
537:       PetscObjectGetComm((PetscObject)vfull,&comm); 
538:       VecGetOwnershipRange(vreduced,&rlow,&rhigh); 
539:       ISCreateStride(comm,rhigh-rlow,rlow,1,&ident); 
540:       VecScatterCreate(vreduced,ident,vfull,is,&scatter); 
541:       VecScatterBegin(scatter,vreduced,vfull,ADD_VALUES,SCATTER_FORWARD); 
542:       VecScatterEnd(scatter,vreduced,vfull,ADD_VALUES,SCATTER_FORWARD); 
543:       VecScatterDestroy(&scatter); 
544:       ISDestroy(&ident); 
545:     }
546:     
547:     return(0);
548: }


553: /*@
554:    ISCreateComplement - Creates the complement of the the index set

556:    Collective on IS

558:    Input Parameter:
559: +  S -  a PETSc IS
560: -  V - the reference vector space

562:    Output Parameter:
563: .  T -  the complement of S


566: .seealso ISCreateGeneral()

568:    Level: advanced
569: @*/
570: PetscErrorCode ISCreateComplement(IS S, Vec V, IS *T){
572:   PetscInt i,nis,nloc,high,low,n=0;
573:   const PetscInt *s;
574:   PetscInt *tt,*ss;
575:   MPI_Comm comm;


581:   VecGetOwnershipRange(V,&low,&high); 
582:   VecGetLocalSize(V,&nloc); 
583:   ISGetLocalSize(S,&nis); 
584:   ISGetIndices(S, &s); 
585:   PetscMalloc( nloc*sizeof(PetscInt),&tt ); 
586:   PetscMalloc( nloc*sizeof(PetscInt),&ss ); 

588:   for (i=low; i<high; i++){ tt[i-low]=i; }

590:   for (i=0; i<nis; i++){ tt[s[i]-low] = -2; }
591:   
592:   for (i=0; i<nloc; i++){
593:     if (tt[i]>-1){ ss[n]=tt[i]; n++; }
594:   }

596:   ISRestoreIndices(S, &s); 
597:   
598:   PetscObjectGetComm((PetscObject)S,&comm);
599:   ISCreateGeneral(comm,n,ss,PETSC_COPY_VALUES,T);
600:   
601:   if (tt) {
602:     PetscFree(tt); 
603:   }
604:   if (ss) {
605:     PetscFree(ss); 
606:   }

608:   return(0);
609: }

613: /*@
614:    VecISSetToConstant - Sets the elements of a vector, specified by an index set, to a constant

616:    Input Parameter:
617: +  S -  a PETSc IS
618: .  c - the constant
619: -  V - a Vec

621: .seealso VecSet()

623:    Level: advanced
624: @*/
625: PetscErrorCode VecISSetToConstant(IS S, PetscReal c, Vec V){
627:   PetscInt nloc,low,high,i;
628:   const PetscInt *s;
629:   PetscReal *v;

636:   VecGetOwnershipRange(V, &low, &high); 
637:   ISGetLocalSize(S,&nloc);

639:   ISGetIndices(S, &s); 
640:   VecGetArray(V,&v); 
641:   for (i=0; i<nloc; i++){
642:     v[s[i]-low] = c;
643:   }
644:   
645:   ISRestoreIndices(S, &s); 
646:   VecRestoreArray(V,&v); 

648:   return(0);

650: }

654: /*@ 
655:   MatGetSubMat - Gets a submatrix using the IS

657:   Input Parameters:
658: + M - the full matrix (n x n)
659: . is - the index set for the submatrix (both row and column index sets need to be the same)
660: . v1 - work vector of dimension n, needed for TAO_SUBSET_MASK option
661: - subset_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK,
662:   TAO_SUBSET_MATRIXFREE)

664:   Output Parameters:
665: . Msub - the submatrix
666: @*/
667: PetscErrorCode MatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
668: {
670:   IS iscomp;
671:   PetscBool flg;
675:   if (*Msub) {
676:     MatDestroy(Msub); 
677:   }
678:   switch (subset_type) {
679:     case TAO_SUBSET_SUBVEC:
680:       MatGetSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub); 
681:       break;

683:     case TAO_SUBSET_MASK:
684:       /* Get Reduced Hessian 
685:          Msub[i,j] = M[i,j] if i,j in Free_Local or i==j
686:          Msub[i,j] = 0      if i!=j and i or j not in Free_Local
687:       */
688:       PetscOptionsBool("-different_submatrix","use separate hessian matrix when computing submatrices","TaoSubsetType",PETSC_FALSE,&flg,PETSC_NULL);
689:       if (flg == PETSC_TRUE) {
690:         MatDuplicate(M, MAT_COPY_VALUES, Msub); 
691:       } else {
692:         /* Act on hessian directly (default) */
693:         PetscObjectReference((PetscObject)M); 
694:         *Msub = M;
695:       }
696:       /* Save the diagonal to temporary vector */
697:       MatGetDiagonal(*Msub,v1); 
698:     
699:       /* Zero out rows and columns */
700:       ISCreateComplement(is,v1,&iscomp); 

702:       /* Use v1 instead of 0 here because of PETSc bug */
703:       MatZeroRowsColumnsIS(*Msub,iscomp,1.0,v1,v1); 

705:       ISDestroy(&iscomp); 
706:       break;
707:     case TAO_SUBSET_MATRIXFREE:
708:       ISCreateComplement(is,v1,&iscomp); 
709:       MatCreateSubMatrixFree(M,iscomp,iscomp,Msub); 
710:       ISDestroy(&iscomp); 
711:       break;
712:   }      
713:   return(0);
714: }