Actual source code: isutil.c
petsc-3.9.4 2018-09-11
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: }