Actual source code: normmh.c
petsc-3.7.5 2017-01-01
2: #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
4: typedef struct {
5: Mat A;
6: Vec w,left,right,leftwork,rightwork;
7: PetscScalar scale;
8: } Mat_Normal;
12: PetscErrorCode MatScaleHermitian_Normal(Mat inA,PetscScalar scale)
13: {
14: Mat_Normal *a = (Mat_Normal*)inA->data;
17: a->scale *= scale;
18: return(0);
19: }
23: PetscErrorCode MatDiagonalScaleHermitian_Normal(Mat inA,Vec left,Vec right)
24: {
25: Mat_Normal *a = (Mat_Normal*)inA->data;
29: if (left) {
30: if (!a->left) {
31: VecDuplicate(left,&a->left);
32: VecCopy(left,a->left);
33: } else {
34: VecPointwiseMult(a->left,left,a->left);
35: }
36: }
37: if (right) {
38: if (!a->right) {
39: VecDuplicate(right,&a->right);
40: VecCopy(right,a->right);
41: } else {
42: VecPointwiseMult(a->right,right,a->right);
43: }
44: }
45: return(0);
46: }
50: PetscErrorCode MatMultHermitian_Normal(Mat N,Vec x,Vec y)
51: {
52: Mat_Normal *Na = (Mat_Normal*)N->data;
54: Vec in;
57: in = x;
58: if (Na->right) {
59: if (!Na->rightwork) {
60: VecDuplicate(Na->right,&Na->rightwork);
61: }
62: VecPointwiseMult(Na->rightwork,Na->right,in);
63: in = Na->rightwork;
64: }
65: MatMult(Na->A,in,Na->w);
66: MatMultHermitianTranspose(Na->A,Na->w,y);
67: if (Na->left) {
68: VecPointwiseMult(y,Na->left,y);
69: }
70: VecScale(y,Na->scale);
71: return(0);
72: }
76: PetscErrorCode MatMultHermitianAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
77: {
78: Mat_Normal *Na = (Mat_Normal*)N->data;
80: Vec in;
83: in = v1;
84: if (Na->right) {
85: if (!Na->rightwork) {
86: VecDuplicate(Na->right,&Na->rightwork);
87: }
88: VecPointwiseMult(Na->rightwork,Na->right,in);
89: in = Na->rightwork;
90: }
91: MatMult(Na->A,in,Na->w);
92: VecScale(Na->w,Na->scale);
93: if (Na->left) {
94: MatMultHermitianTranspose(Na->A,Na->w,v3);
95: VecPointwiseMult(v3,Na->left,v3);
96: VecAXPY(v3,1.0,v2);
97: } else {
98: MatMultHermitianTransposeAdd(Na->A,Na->w,v2,v3);
99: }
100: return(0);
101: }
105: PetscErrorCode MatMultHermitianTranspose_Normal(Mat N,Vec x,Vec y)
106: {
107: Mat_Normal *Na = (Mat_Normal*)N->data;
109: Vec in;
112: in = x;
113: if (Na->left) {
114: if (!Na->leftwork) {
115: VecDuplicate(Na->left,&Na->leftwork);
116: }
117: VecPointwiseMult(Na->leftwork,Na->left,in);
118: in = Na->leftwork;
119: }
120: MatMult(Na->A,in,Na->w);
121: MatMultHermitianTranspose(Na->A,Na->w,y);
122: if (Na->right) {
123: VecPointwiseMult(y,Na->right,y);
124: }
125: VecScale(y,Na->scale);
126: return(0);
127: }
131: PetscErrorCode MatMultHermitianTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
132: {
133: Mat_Normal *Na = (Mat_Normal*)N->data;
135: Vec in;
138: in = v1;
139: if (Na->left) {
140: if (!Na->leftwork) {
141: VecDuplicate(Na->left,&Na->leftwork);
142: }
143: VecPointwiseMult(Na->leftwork,Na->left,in);
144: in = Na->leftwork;
145: }
146: MatMult(Na->A,in,Na->w);
147: VecScale(Na->w,Na->scale);
148: if (Na->right) {
149: MatMultHermitianTranspose(Na->A,Na->w,v3);
150: VecPointwiseMult(v3,Na->right,v3);
151: VecAXPY(v3,1.0,v2);
152: } else {
153: MatMultHermitianTransposeAdd(Na->A,Na->w,v2,v3);
154: }
155: return(0);
156: }
160: PetscErrorCode MatDestroyHermitian_Normal(Mat N)
161: {
162: Mat_Normal *Na = (Mat_Normal*)N->data;
166: MatDestroy(&Na->A);
167: VecDestroy(&Na->w);
168: VecDestroy(&Na->left);
169: VecDestroy(&Na->right);
170: VecDestroy(&Na->leftwork);
171: VecDestroy(&Na->rightwork);
172: PetscFree(N->data);
173: return(0);
174: }
176: /*
177: Slow, nonscalable version
178: */
181: PetscErrorCode MatGetDiagonalHermitian_Normal(Mat N,Vec v)
182: {
183: Mat_Normal *Na = (Mat_Normal*)N->data;
184: Mat A = Na->A;
185: PetscErrorCode ierr;
186: PetscInt i,j,rstart,rend,nnz;
187: const PetscInt *cols;
188: PetscScalar *diag,*work,*values;
189: const PetscScalar *mvalues;
192: PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);
193: PetscMemzero(work,A->cmap->N*sizeof(PetscScalar));
194: MatGetOwnershipRange(A,&rstart,&rend);
195: for (i=rstart; i<rend; i++) {
196: MatGetRow(A,i,&nnz,&cols,&mvalues);
197: for (j=0; j<nnz; j++) {
198: work[cols[j]] += mvalues[j]*PetscConj(mvalues[j]);
199: }
200: MatRestoreRow(A,i,&nnz,&cols,&mvalues);
201: }
202: MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));
203: rstart = N->cmap->rstart;
204: rend = N->cmap->rend;
205: VecGetArray(v,&values);
206: PetscMemcpy(values,diag+rstart,(rend-rstart)*sizeof(PetscScalar));
207: VecRestoreArray(v,&values);
208: PetscFree2(diag,work);
209: VecScale(v,Na->scale);
210: return(0);
211: }
215: /*@
216: MatCreateNormalHermitian - Creates a new matrix object that behaves like (A*)'*A.
218: Collective on Mat
220: Input Parameter:
221: . A - the (possibly rectangular complex) matrix
223: Output Parameter:
224: . N - the matrix that represents (A*)'*A
226: Level: intermediate
228: Notes: The product (A*)'*A is NOT actually formed! Rather the new matrix
229: object performs the matrix-vector product by first multiplying by
230: A and then (A*)'
231: @*/
232: PetscErrorCode MatCreateNormalHermitian(Mat A,Mat *N)
233: {
235: PetscInt m,n;
236: Mat_Normal *Na;
239: MatGetLocalSize(A,&m,&n);
240: MatCreate(PetscObjectComm((PetscObject)A),N);
241: MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);
242: PetscObjectChangeTypeName((PetscObject)*N,MATNORMALHERMITIAN);
244: PetscNewLog(*N,&Na);
245: (*N)->data = (void*) Na;
246: PetscObjectReference((PetscObject)A);
247: Na->A = A;
248: Na->scale = 1.0;
250: VecCreateMPI(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,&Na->w);
252: (*N)->ops->destroy = MatDestroyHermitian_Normal;
253: (*N)->ops->mult = MatMultHermitian_Normal;
254: (*N)->ops->multtranspose = MatMultHermitianTranspose_Normal;
255: (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_Normal;
256: (*N)->ops->multadd = MatMultHermitianAdd_Normal;
257: (*N)->ops->getdiagonal = MatGetDiagonalHermitian_Normal;
258: (*N)->ops->scale = MatScaleHermitian_Normal;
259: (*N)->ops->diagonalscale = MatDiagonalScaleHermitian_Normal;
260: (*N)->assembled = PETSC_TRUE;
261: (*N)->cmap->N = A->cmap->N;
262: (*N)->rmap->N = A->cmap->N;
263: (*N)->cmap->n = A->cmap->n;
264: (*N)->rmap->n = A->cmap->n;
265: return(0);
266: }