1:  #include <petsctao.h>
  2:  #include <petsc/private/taoimpl.h>
  3:  #include <../src/tao/matrix/submatfree.h>

  5: /*@C
  6:   TaoVecGetSubVec - Gets a subvector using the IS

  8:   Input Parameters:
  9: + vfull - the full matrix
 10: . is - the index set for the subvector
 11: . reduced_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK,  TAO_SUBSET_MATRIXFREE)
 12: - maskvalue - the value to set the unused vector elements to (for TAO_SUBSET_MASK or TAO_SUBSET_MATRIXFREE)

 14:   Output Parameters:
 15: . vreduced - the subvector

 17:   Notes:
 18:   maskvalue should usually be 0.0, unless a pointwise divide will be used.

 20: @*/
 21: PetscErrorCode TaoVecGetSubVec(Vec vfull, IS is, TaoSubsetType reduced_type, PetscReal maskvalue, Vec *vreduced)
 22: {
 24:   PetscInt       nfull,nreduced,nreduced_local,rlow,rhigh,flow,fhigh;
 25:   PetscInt       i,nlocal;
 26:   PetscReal      *fv,*rv;
 27:   const PetscInt *s;
 28:   IS             ident;
 29:   VecType        vtype;
 30:   VecScatter     scatter;
 31:   MPI_Comm       comm;


 37:   VecGetSize(vfull, &nfull);
 38:   ISGetSize(is, &nreduced);

 40:   if (nreduced == nfull) {
 41:     VecDestroy(vreduced);
 42:     VecDuplicate(vfull,vreduced);
 43:     VecCopy(vfull,*vreduced);
 44:   } else {
 45:     switch (reduced_type) {
 46:     case TAO_SUBSET_SUBVEC:
 47:       VecGetType(vfull,&vtype);
 48:       VecGetOwnershipRange(vfull,&flow,&fhigh);
 49:       ISGetLocalSize(is,&nreduced_local);
 50:       PetscObjectGetComm((PetscObject)vfull,&comm);
 51:       if (*vreduced) {
 52:         VecDestroy(vreduced);
 53:       }
 54:       VecCreate(comm,vreduced);
 55:       VecSetType(*vreduced,vtype);

 57:       VecSetSizes(*vreduced,nreduced_local,nreduced);
 58:       VecGetOwnershipRange(*vreduced,&rlow,&rhigh);
 59:       ISCreateStride(comm,nreduced_local,rlow,1,&ident);
 60:       VecScatterCreate(vfull,is,*vreduced,ident,&scatter);
 61:       VecScatterBegin(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);
 62:       VecScatterEnd(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);
 63:       VecScatterDestroy(&scatter);
 64:       ISDestroy(&ident);
 65:       break;

 67:     case TAO_SUBSET_MASK:
 68:     case TAO_SUBSET_MATRIXFREE:
 69:       /* vr[i] = vf[i]   if i in is
 70:        vr[i] = 0       otherwise */
 71:       if (!*vreduced) {
 72:         VecDuplicate(vfull,vreduced);
 73:       }

 75:       VecSet(*vreduced,maskvalue);
 76:       ISGetLocalSize(is,&nlocal);
 77:       VecGetOwnershipRange(vfull,&flow,&fhigh);
 78:       VecGetArray(vfull,&fv);
 79:       VecGetArray(*vreduced,&rv);
 80:       ISGetIndices(is,&s);
 81:       if (nlocal > (fhigh-flow)) SETERRQ2(PETSC_COMM_WORLD,1,"IS local size %d > Vec local size %d",nlocal,fhigh-flow);
 82:       for (i=0;i<nlocal;i++) {
 83:         rv[s[i]-flow] = fv[s[i]-flow];
 84:       }
 85:       ISRestoreIndices(is,&s);
 86:       VecRestoreArray(vfull,&fv);
 87:       VecRestoreArray(*vreduced,&rv);
 88:       break;
 89:     }
 90:   }
 91:   return(0);
 92: }

 94: /*@C
 95:   TaoMatGetSubMat - Gets a submatrix using the IS

 97:   Input Parameters:
 98: + M - the full matrix (n x n)
 99: . is - the index set for the submatrix (both row and column index sets need to be the same)
100: . v1 - work vector of dimension n, needed for TAO_SUBSET_MASK option
101: - subset_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK,
102:   TAO_SUBSET_MATRIXFREE)

104:   Output Parameters:
105: . Msub - the submatrix
106: @*/
107: PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
108: {
110:   IS             iscomp;
111:   PetscBool      flg = PETSC_FALSE;

116:   MatDestroy(Msub);
117:   switch (subset_type) {
118:   case TAO_SUBSET_SUBVEC:
119:     MatCreateSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);
120:     break;

122:   case TAO_SUBSET_MASK:
123:     /* Get Reduced Hessian
124:      Msub[i,j] = M[i,j] if i,j in Free_Local or i==j
125:      Msub[i,j] = 0      if i!=j and i or j not in Free_Local
126:      */
127:     PetscOptionsBegin(PetscObjectComm((PetscObject)M),NULL,NULL,NULL);
128:     PetscOptionsBool("-different_submatrix","use separate hessian matrix when computing submatrices","TaoSubsetType",flg,&flg,NULL);
129:     PetscOptionsEnd();
130:     if (flg) {
131:       MatDuplicate(M, MAT_COPY_VALUES, Msub);
132:     } else {
133:       /* Act on hessian directly (default) */
134:       PetscObjectReference((PetscObject)M);
135:       *Msub = M;
136:     }
137:     /* Save the diagonal to temporary vector */
138:     MatGetDiagonal(*Msub,v1);

140:     /* Zero out rows and columns */
141:     ISComplementVec(is,v1,&iscomp);

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

146:     ISDestroy(&iscomp);
147:     break;
148:   case TAO_SUBSET_MATRIXFREE:
149:     ISComplementVec(is,v1,&iscomp);
150:     MatCreateSubMatrixFree(M,iscomp,iscomp,Msub);
151:     ISDestroy(&iscomp);
152:     break;
153:   }
154:   return(0);
155: }