2: #include <petsc/private/matimpl.h>

  4: typedef struct {
  5:   Mat         A;
  6:   Vec         w,left,right,leftwork,rightwork;
  7:   PetscScalar scale;
  8: } Mat_Normal;

 10: PetscErrorCode MatScaleHermitian_Normal(Mat inA,PetscScalar scale)
 11: {
 12:   Mat_Normal *a = (Mat_Normal*)inA->data;

 14:   a->scale *= scale;
 15:   return 0;
 16: }

 18: PetscErrorCode MatDiagonalScaleHermitian_Normal(Mat inA,Vec left,Vec right)
 19: {
 20:   Mat_Normal     *a = (Mat_Normal*)inA->data;

 22:   if (left) {
 23:     if (!a->left) {
 24:       VecDuplicate(left,&a->left);
 25:       VecCopy(left,a->left);
 26:     } else {
 27:       VecPointwiseMult(a->left,left,a->left);
 28:     }
 29:   }
 30:   if (right) {
 31:     if (!a->right) {
 32:       VecDuplicate(right,&a->right);
 33:       VecCopy(right,a->right);
 34:     } else {
 35:       VecPointwiseMult(a->right,right,a->right);
 36:     }
 37:   }
 38:   return 0;
 39: }

 41: PetscErrorCode MatMultHermitian_Normal(Mat N,Vec x,Vec y)
 42: {
 43:   Mat_Normal     *Na = (Mat_Normal*)N->data;
 44:   Vec            in;

 46:   in = x;
 47:   if (Na->right) {
 48:     if (!Na->rightwork) {
 49:       VecDuplicate(Na->right,&Na->rightwork);
 50:     }
 51:     VecPointwiseMult(Na->rightwork,Na->right,in);
 52:     in   = Na->rightwork;
 53:   }
 54:   MatMult(Na->A,in,Na->w);
 55:   MatMultHermitianTranspose(Na->A,Na->w,y);
 56:   if (Na->left) {
 57:     VecPointwiseMult(y,Na->left,y);
 58:   }
 59:   VecScale(y,Na->scale);
 60:   return 0;
 61: }

 63: PetscErrorCode MatMultHermitianAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
 64: {
 65:   Mat_Normal     *Na = (Mat_Normal*)N->data;
 66:   Vec            in;

 68:   in = v1;
 69:   if (Na->right) {
 70:     if (!Na->rightwork) {
 71:       VecDuplicate(Na->right,&Na->rightwork);
 72:     }
 73:     VecPointwiseMult(Na->rightwork,Na->right,in);
 74:     in   = Na->rightwork;
 75:   }
 76:   MatMult(Na->A,in,Na->w);
 77:   VecScale(Na->w,Na->scale);
 78:   if (Na->left) {
 79:     MatMultHermitianTranspose(Na->A,Na->w,v3);
 80:     VecPointwiseMult(v3,Na->left,v3);
 81:     VecAXPY(v3,1.0,v2);
 82:   } else {
 83:     MatMultHermitianTransposeAdd(Na->A,Na->w,v2,v3);
 84:   }
 85:   return 0;
 86: }

 88: PetscErrorCode MatMultHermitianTranspose_Normal(Mat N,Vec x,Vec y)
 89: {
 90:   Mat_Normal     *Na = (Mat_Normal*)N->data;
 91:   Vec            in;

 93:   in = x;
 94:   if (Na->left) {
 95:     if (!Na->leftwork) {
 96:       VecDuplicate(Na->left,&Na->leftwork);
 97:     }
 98:     VecPointwiseMult(Na->leftwork,Na->left,in);
 99:     in   = Na->leftwork;
100:   }
101:   MatMult(Na->A,in,Na->w);
102:   MatMultHermitianTranspose(Na->A,Na->w,y);
103:   if (Na->right) {
104:     VecPointwiseMult(y,Na->right,y);
105:   }
106:   VecScale(y,Na->scale);
107:   return 0;
108: }

110: PetscErrorCode MatMultHermitianTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
111: {
112:   Mat_Normal     *Na = (Mat_Normal*)N->data;
113:   Vec            in;

115:   in = v1;
116:   if (Na->left) {
117:     if (!Na->leftwork) {
118:       VecDuplicate(Na->left,&Na->leftwork);
119:     }
120:     VecPointwiseMult(Na->leftwork,Na->left,in);
121:     in   = Na->leftwork;
122:   }
123:   MatMult(Na->A,in,Na->w);
124:   VecScale(Na->w,Na->scale);
125:   if (Na->right) {
126:     MatMultHermitianTranspose(Na->A,Na->w,v3);
127:     VecPointwiseMult(v3,Na->right,v3);
128:     VecAXPY(v3,1.0,v2);
129:   } else {
130:     MatMultHermitianTransposeAdd(Na->A,Na->w,v2,v3);
131:   }
132:   return 0;
133: }

135: PetscErrorCode MatDestroyHermitian_Normal(Mat N)
136: {
137:   Mat_Normal     *Na = (Mat_Normal*)N->data;

139:   MatDestroy(&Na->A);
140:   VecDestroy(&Na->w);
141:   VecDestroy(&Na->left);
142:   VecDestroy(&Na->right);
143:   VecDestroy(&Na->leftwork);
144:   VecDestroy(&Na->rightwork);
145:   PetscFree(N->data);
146:   PetscObjectComposeFunction((PetscObject)N,"MatNormalGetMatHermitian_C",NULL);
147:   return 0;
148: }

150: /*
151:       Slow, nonscalable version
152: */
153: PetscErrorCode MatGetDiagonalHermitian_Normal(Mat N,Vec v)
154: {
155:   Mat_Normal        *Na = (Mat_Normal*)N->data;
156:   Mat               A   = Na->A;
157:   PetscInt          i,j,rstart,rend,nnz;
158:   const PetscInt    *cols;
159:   PetscScalar       *diag,*work,*values;
160:   const PetscScalar *mvalues;

162:   PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);
163:   PetscArrayzero(work,A->cmap->N);
164:   MatGetOwnershipRange(A,&rstart,&rend);
165:   for (i=rstart; i<rend; i++) {
166:     MatGetRow(A,i,&nnz,&cols,&mvalues);
167:     for (j=0; j<nnz; j++) {
168:       work[cols[j]] += mvalues[j]*PetscConj(mvalues[j]);
169:     }
170:     MatRestoreRow(A,i,&nnz,&cols,&mvalues);
171:   }
172:   MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));
173:   rstart = N->cmap->rstart;
174:   rend   = N->cmap->rend;
175:   VecGetArray(v,&values);
176:   PetscArraycpy(values,diag+rstart,rend-rstart);
177:   VecRestoreArray(v,&values);
178:   PetscFree2(diag,work);
179:   VecScale(v,Na->scale);
180:   return 0;
181: }

183: PetscErrorCode MatNormalGetMatHermitian_Normal(Mat A,Mat *M)
184: {
185:   Mat_Normal *Aa = (Mat_Normal*)A->data;

187:   *M = Aa->A;
188:   return 0;
189: }

191: /*@
192:       MatNormalHermitianGetMat - Gets the Mat object stored inside a MATNORMALHERMITIAN

194:    Logically collective on Mat

196:    Input Parameter:
197: .   A  - the MATNORMALHERMITIAN matrix

199:    Output Parameter:
200: .   M - the matrix object stored inside A

202:    Level: intermediate

204: .seealso: MatCreateNormalHermitian()

206: @*/
207: PetscErrorCode MatNormalHermitianGetMat(Mat A,Mat *M)
208: {
212:   PetscUseMethod(A,"MatNormalGetMatHermitian_C",(Mat,Mat*),(A,M));
213:   return 0;
214: }

216: /*@
217:       MatCreateNormalHermitian - Creates a new matrix object that behaves like (A*)'*A.

219:    Collective on Mat

221:    Input Parameter:
222: .   A  - the (possibly rectangular complex) matrix

224:    Output Parameter:
225: .   N - the matrix that represents (A*)'*A

227:    Level: intermediate

229:    Notes:
230:     The product (A*)'*A is NOT actually formed! Rather the new matrix
231:           object performs the matrix-vector product by first multiplying by
232:           A and then (A*)'
233: @*/
234: PetscErrorCode  MatCreateNormalHermitian(Mat A,Mat *N)
235: {
236:   PetscInt       m,n;
237:   Mat_Normal     *Na;
238:   VecType        vtype;

240:   MatGetLocalSize(A,&m,&n);
241:   MatCreate(PetscObjectComm((PetscObject)A),N);
242:   MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);
243:   PetscObjectChangeTypeName((PetscObject)*N,MATNORMALHERMITIAN);
244:   PetscLayoutReference(A->cmap,&(*N)->rmap);
245:   PetscLayoutReference(A->cmap,&(*N)->cmap);

247:   PetscNewLog(*N,&Na);
248:   (*N)->data = (void*) Na;
249:   PetscObjectReference((PetscObject)A);
250:   Na->A      = A;
251:   Na->scale  = 1.0;

253:   MatCreateVecs(A,NULL,&Na->w);

255:   (*N)->ops->destroy          = MatDestroyHermitian_Normal;
256:   (*N)->ops->mult             = MatMultHermitian_Normal;
257:   (*N)->ops->multtranspose    = MatMultHermitianTranspose_Normal;
258:   (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_Normal;
259:   (*N)->ops->multadd          = MatMultHermitianAdd_Normal;
260:   (*N)->ops->getdiagonal      = MatGetDiagonalHermitian_Normal;
261:   (*N)->ops->scale            = MatScaleHermitian_Normal;
262:   (*N)->ops->diagonalscale    = MatDiagonalScaleHermitian_Normal;
263:   (*N)->assembled             = PETSC_TRUE;
264:   (*N)->preallocated          = PETSC_TRUE;

266:   PetscObjectComposeFunction((PetscObject)(*N),"MatNormalGetMatHermitian_C",MatNormalGetMatHermitian_Normal);
267:   MatSetOption(*N,MAT_HERMITIAN,PETSC_TRUE);
268:   MatGetVecType(A,&vtype);
269:   MatSetVecType(*N,vtype);
270: #if defined(PETSC_HAVE_DEVICE)
271:   MatBindToCPU(*N,A->boundtocpu);
272: #endif
273:   return 0;
274: }