Actual source code: normmh.c
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;
15: a->scale *= scale;
16: return(0);
17: }
19: PetscErrorCode MatDiagonalScaleHermitian_Normal(Mat inA,Vec left,Vec right)
20: {
21: Mat_Normal *a = (Mat_Normal*)inA->data;
25: if (left) {
26: if (!a->left) {
27: VecDuplicate(left,&a->left);
28: VecCopy(left,a->left);
29: } else {
30: VecPointwiseMult(a->left,left,a->left);
31: }
32: }
33: if (right) {
34: if (!a->right) {
35: VecDuplicate(right,&a->right);
36: VecCopy(right,a->right);
37: } else {
38: VecPointwiseMult(a->right,right,a->right);
39: }
40: }
41: return(0);
42: }
44: PetscErrorCode MatMultHermitian_Normal(Mat N,Vec x,Vec y)
45: {
46: Mat_Normal *Na = (Mat_Normal*)N->data;
48: Vec in;
51: in = x;
52: if (Na->right) {
53: if (!Na->rightwork) {
54: VecDuplicate(Na->right,&Na->rightwork);
55: }
56: VecPointwiseMult(Na->rightwork,Na->right,in);
57: in = Na->rightwork;
58: }
59: MatMult(Na->A,in,Na->w);
60: MatMultHermitianTranspose(Na->A,Na->w,y);
61: if (Na->left) {
62: VecPointwiseMult(y,Na->left,y);
63: }
64: VecScale(y,Na->scale);
65: return(0);
66: }
68: PetscErrorCode MatMultHermitianAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
69: {
70: Mat_Normal *Na = (Mat_Normal*)N->data;
72: Vec in;
75: in = v1;
76: if (Na->right) {
77: if (!Na->rightwork) {
78: VecDuplicate(Na->right,&Na->rightwork);
79: }
80: VecPointwiseMult(Na->rightwork,Na->right,in);
81: in = Na->rightwork;
82: }
83: MatMult(Na->A,in,Na->w);
84: VecScale(Na->w,Na->scale);
85: if (Na->left) {
86: MatMultHermitianTranspose(Na->A,Na->w,v3);
87: VecPointwiseMult(v3,Na->left,v3);
88: VecAXPY(v3,1.0,v2);
89: } else {
90: MatMultHermitianTransposeAdd(Na->A,Na->w,v2,v3);
91: }
92: return(0);
93: }
95: PetscErrorCode MatMultHermitianTranspose_Normal(Mat N,Vec x,Vec y)
96: {
97: Mat_Normal *Na = (Mat_Normal*)N->data;
99: Vec in;
102: in = x;
103: if (Na->left) {
104: if (!Na->leftwork) {
105: VecDuplicate(Na->left,&Na->leftwork);
106: }
107: VecPointwiseMult(Na->leftwork,Na->left,in);
108: in = Na->leftwork;
109: }
110: MatMult(Na->A,in,Na->w);
111: MatMultHermitianTranspose(Na->A,Na->w,y);
112: if (Na->right) {
113: VecPointwiseMult(y,Na->right,y);
114: }
115: VecScale(y,Na->scale);
116: return(0);
117: }
119: PetscErrorCode MatMultHermitianTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
120: {
121: Mat_Normal *Na = (Mat_Normal*)N->data;
123: Vec in;
126: in = v1;
127: if (Na->left) {
128: if (!Na->leftwork) {
129: VecDuplicate(Na->left,&Na->leftwork);
130: }
131: VecPointwiseMult(Na->leftwork,Na->left,in);
132: in = Na->leftwork;
133: }
134: MatMult(Na->A,in,Na->w);
135: VecScale(Na->w,Na->scale);
136: if (Na->right) {
137: MatMultHermitianTranspose(Na->A,Na->w,v3);
138: VecPointwiseMult(v3,Na->right,v3);
139: VecAXPY(v3,1.0,v2);
140: } else {
141: MatMultHermitianTransposeAdd(Na->A,Na->w,v2,v3);
142: }
143: return(0);
144: }
146: PetscErrorCode MatDestroyHermitian_Normal(Mat N)
147: {
148: Mat_Normal *Na = (Mat_Normal*)N->data;
152: MatDestroy(&Na->A);
153: VecDestroy(&Na->w);
154: VecDestroy(&Na->left);
155: VecDestroy(&Na->right);
156: VecDestroy(&Na->leftwork);
157: VecDestroy(&Na->rightwork);
158: PetscFree(N->data);
159: PetscObjectComposeFunction((PetscObject)N,"MatNormalGetMatHermitian_C",NULL);
160: return(0);
161: }
163: /*
164: Slow, nonscalable version
165: */
166: PetscErrorCode MatGetDiagonalHermitian_Normal(Mat N,Vec v)
167: {
168: Mat_Normal *Na = (Mat_Normal*)N->data;
169: Mat A = Na->A;
170: PetscErrorCode ierr;
171: PetscInt i,j,rstart,rend,nnz;
172: const PetscInt *cols;
173: PetscScalar *diag,*work,*values;
174: const PetscScalar *mvalues;
177: PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);
178: PetscArrayzero(work,A->cmap->N);
179: MatGetOwnershipRange(A,&rstart,&rend);
180: for (i=rstart; i<rend; i++) {
181: MatGetRow(A,i,&nnz,&cols,&mvalues);
182: for (j=0; j<nnz; j++) {
183: work[cols[j]] += mvalues[j]*PetscConj(mvalues[j]);
184: }
185: MatRestoreRow(A,i,&nnz,&cols,&mvalues);
186: }
187: MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));
188: rstart = N->cmap->rstart;
189: rend = N->cmap->rend;
190: VecGetArray(v,&values);
191: PetscArraycpy(values,diag+rstart,rend-rstart);
192: VecRestoreArray(v,&values);
193: PetscFree2(diag,work);
194: VecScale(v,Na->scale);
195: return(0);
196: }
198: PetscErrorCode MatNormalGetMatHermitian_Normal(Mat A,Mat *M)
199: {
200: Mat_Normal *Aa = (Mat_Normal*)A->data;
203: *M = Aa->A;
204: return(0);
205: }
207: /*@
208: MatNormalHermitianGetMat - Gets the Mat object stored inside a MATNORMALHERMITIAN
210: Logically collective on Mat
212: Input Parameter:
213: . A - the MATNORMALHERMITIAN matrix
215: Output Parameter:
216: . M - the matrix object stored inside A
218: Level: intermediate
220: .seealso: MatCreateNormalHermitian()
222: @*/
223: PetscErrorCode MatNormalHermitianGetMat(Mat A,Mat *M)
224: {
231: PetscUseMethod(A,"MatNormalGetMatHermitian_C",(Mat,Mat*),(A,M));
232: return(0);
233: }
235: /*@
236: MatCreateNormalHermitian - Creates a new matrix object that behaves like (A*)'*A.
238: Collective on Mat
240: Input Parameter:
241: . A - the (possibly rectangular complex) matrix
243: Output Parameter:
244: . N - the matrix that represents (A*)'*A
246: Level: intermediate
248: Notes:
249: The product (A*)'*A is NOT actually formed! Rather the new matrix
250: object performs the matrix-vector product by first multiplying by
251: A and then (A*)'
252: @*/
253: PetscErrorCode MatCreateNormalHermitian(Mat A,Mat *N)
254: {
256: PetscInt m,n;
257: Mat_Normal *Na;
258: VecType vtype;
261: MatGetLocalSize(A,&m,&n);
262: MatCreate(PetscObjectComm((PetscObject)A),N);
263: MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);
264: PetscObjectChangeTypeName((PetscObject)*N,MATNORMALHERMITIAN);
265: PetscLayoutReference(A->cmap,&(*N)->rmap);
266: PetscLayoutReference(A->cmap,&(*N)->cmap);
268: PetscNewLog(*N,&Na);
269: (*N)->data = (void*) Na;
270: PetscObjectReference((PetscObject)A);
271: Na->A = A;
272: Na->scale = 1.0;
274: MatCreateVecs(A,NULL,&Na->w);
276: (*N)->ops->destroy = MatDestroyHermitian_Normal;
277: (*N)->ops->mult = MatMultHermitian_Normal;
278: (*N)->ops->multtranspose = MatMultHermitianTranspose_Normal;
279: (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_Normal;
280: (*N)->ops->multadd = MatMultHermitianAdd_Normal;
281: (*N)->ops->getdiagonal = MatGetDiagonalHermitian_Normal;
282: (*N)->ops->scale = MatScaleHermitian_Normal;
283: (*N)->ops->diagonalscale = MatDiagonalScaleHermitian_Normal;
284: (*N)->assembled = PETSC_TRUE;
285: (*N)->preallocated = PETSC_TRUE;
287: PetscObjectComposeFunction((PetscObject)(*N),"MatNormalGetMatHermitian_C",MatNormalGetMatHermitian_Normal);
288: MatSetOption(*N,MAT_HERMITIAN,PETSC_TRUE);
289: MatGetVecType(A,&vtype);
290: MatSetVecType(*N,vtype);
291: #if defined(PETSC_HAVE_DEVICE)
292: MatBindToCPU(*N,A->boundtocpu);
293: #endif
294: return(0);
295: }