1
2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
3
4 typedef struct {
5 Mat A;
6 Vec w,left,right,leftwork,rightwork;
7 PetscScalar scale;
8 } Mat_Normal;
9
MatScaleHermitian_Normal(Mat inA,PetscScalar scale)10 PetscErrorCode MatScaleHermitian_Normal(Mat inA,PetscScalar scale)
11 {
12 Mat_Normal *a = (Mat_Normal*)inA->data;
13
14 PetscFunctionBegin;
15 a->scale *= scale;
16 PetscFunctionReturn(0);
17 }
18
MatDiagonalScaleHermitian_Normal(Mat inA,Vec left,Vec right)19 PetscErrorCode MatDiagonalScaleHermitian_Normal(Mat inA,Vec left,Vec right)
20 {
21 Mat_Normal *a = (Mat_Normal*)inA->data;
22 PetscErrorCode ierr;
23
24 PetscFunctionBegin;
25 if (left) {
26 if (!a->left) {
27 ierr = VecDuplicate(left,&a->left);CHKERRQ(ierr);
28 ierr = VecCopy(left,a->left);CHKERRQ(ierr);
29 } else {
30 ierr = VecPointwiseMult(a->left,left,a->left);CHKERRQ(ierr);
31 }
32 }
33 if (right) {
34 if (!a->right) {
35 ierr = VecDuplicate(right,&a->right);CHKERRQ(ierr);
36 ierr = VecCopy(right,a->right);CHKERRQ(ierr);
37 } else {
38 ierr = VecPointwiseMult(a->right,right,a->right);CHKERRQ(ierr);
39 }
40 }
41 PetscFunctionReturn(0);
42 }
43
MatMultHermitian_Normal(Mat N,Vec x,Vec y)44 PetscErrorCode MatMultHermitian_Normal(Mat N,Vec x,Vec y)
45 {
46 Mat_Normal *Na = (Mat_Normal*)N->data;
47 PetscErrorCode ierr;
48 Vec in;
49
50 PetscFunctionBegin;
51 in = x;
52 if (Na->right) {
53 if (!Na->rightwork) {
54 ierr = VecDuplicate(Na->right,&Na->rightwork);CHKERRQ(ierr);
55 }
56 ierr = VecPointwiseMult(Na->rightwork,Na->right,in);CHKERRQ(ierr);
57 in = Na->rightwork;
58 }
59 ierr = MatMult(Na->A,in,Na->w);CHKERRQ(ierr);
60 ierr = MatMultHermitianTranspose(Na->A,Na->w,y);CHKERRQ(ierr);
61 if (Na->left) {
62 ierr = VecPointwiseMult(y,Na->left,y);CHKERRQ(ierr);
63 }
64 ierr = VecScale(y,Na->scale);CHKERRQ(ierr);
65 PetscFunctionReturn(0);
66 }
67
MatMultHermitianAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)68 PetscErrorCode MatMultHermitianAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
69 {
70 Mat_Normal *Na = (Mat_Normal*)N->data;
71 PetscErrorCode ierr;
72 Vec in;
73
74 PetscFunctionBegin;
75 in = v1;
76 if (Na->right) {
77 if (!Na->rightwork) {
78 ierr = VecDuplicate(Na->right,&Na->rightwork);CHKERRQ(ierr);
79 }
80 ierr = VecPointwiseMult(Na->rightwork,Na->right,in);CHKERRQ(ierr);
81 in = Na->rightwork;
82 }
83 ierr = MatMult(Na->A,in,Na->w);CHKERRQ(ierr);
84 ierr = VecScale(Na->w,Na->scale);CHKERRQ(ierr);
85 if (Na->left) {
86 ierr = MatMultHermitianTranspose(Na->A,Na->w,v3);CHKERRQ(ierr);
87 ierr = VecPointwiseMult(v3,Na->left,v3);CHKERRQ(ierr);
88 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr);
89 } else {
90 ierr = MatMultHermitianTransposeAdd(Na->A,Na->w,v2,v3);CHKERRQ(ierr);
91 }
92 PetscFunctionReturn(0);
93 }
94
MatMultHermitianTranspose_Normal(Mat N,Vec x,Vec y)95 PetscErrorCode MatMultHermitianTranspose_Normal(Mat N,Vec x,Vec y)
96 {
97 Mat_Normal *Na = (Mat_Normal*)N->data;
98 PetscErrorCode ierr;
99 Vec in;
100
101 PetscFunctionBegin;
102 in = x;
103 if (Na->left) {
104 if (!Na->leftwork) {
105 ierr = VecDuplicate(Na->left,&Na->leftwork);CHKERRQ(ierr);
106 }
107 ierr = VecPointwiseMult(Na->leftwork,Na->left,in);CHKERRQ(ierr);
108 in = Na->leftwork;
109 }
110 ierr = MatMult(Na->A,in,Na->w);CHKERRQ(ierr);
111 ierr = MatMultHermitianTranspose(Na->A,Na->w,y);CHKERRQ(ierr);
112 if (Na->right) {
113 ierr = VecPointwiseMult(y,Na->right,y);CHKERRQ(ierr);
114 }
115 ierr = VecScale(y,Na->scale);CHKERRQ(ierr);
116 PetscFunctionReturn(0);
117 }
118
MatMultHermitianTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)119 PetscErrorCode MatMultHermitianTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
120 {
121 Mat_Normal *Na = (Mat_Normal*)N->data;
122 PetscErrorCode ierr;
123 Vec in;
124
125 PetscFunctionBegin;
126 in = v1;
127 if (Na->left) {
128 if (!Na->leftwork) {
129 ierr = VecDuplicate(Na->left,&Na->leftwork);CHKERRQ(ierr);
130 }
131 ierr = VecPointwiseMult(Na->leftwork,Na->left,in);CHKERRQ(ierr);
132 in = Na->leftwork;
133 }
134 ierr = MatMult(Na->A,in,Na->w);CHKERRQ(ierr);
135 ierr = VecScale(Na->w,Na->scale);CHKERRQ(ierr);
136 if (Na->right) {
137 ierr = MatMultHermitianTranspose(Na->A,Na->w,v3);CHKERRQ(ierr);
138 ierr = VecPointwiseMult(v3,Na->right,v3);CHKERRQ(ierr);
139 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr);
140 } else {
141 ierr = MatMultHermitianTransposeAdd(Na->A,Na->w,v2,v3);CHKERRQ(ierr);
142 }
143 PetscFunctionReturn(0);
144 }
145
MatDestroyHermitian_Normal(Mat N)146 PetscErrorCode MatDestroyHermitian_Normal(Mat N)
147 {
148 Mat_Normal *Na = (Mat_Normal*)N->data;
149 PetscErrorCode ierr;
150
151 PetscFunctionBegin;
152 ierr = MatDestroy(&Na->A);CHKERRQ(ierr);
153 ierr = VecDestroy(&Na->w);CHKERRQ(ierr);
154 ierr = VecDestroy(&Na->left);CHKERRQ(ierr);
155 ierr = VecDestroy(&Na->right);CHKERRQ(ierr);
156 ierr = VecDestroy(&Na->leftwork);CHKERRQ(ierr);
157 ierr = VecDestroy(&Na->rightwork);CHKERRQ(ierr);
158 ierr = PetscFree(N->data);CHKERRQ(ierr);
159 PetscFunctionReturn(0);
160 }
161
162 /*
163 Slow, nonscalable version
164 */
MatGetDiagonalHermitian_Normal(Mat N,Vec v)165 PetscErrorCode MatGetDiagonalHermitian_Normal(Mat N,Vec v)
166 {
167 Mat_Normal *Na = (Mat_Normal*)N->data;
168 Mat A = Na->A;
169 PetscErrorCode ierr;
170 PetscInt i,j,rstart,rend,nnz;
171 const PetscInt *cols;
172 PetscScalar *diag,*work,*values;
173 const PetscScalar *mvalues;
174
175 PetscFunctionBegin;
176 ierr = PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);CHKERRQ(ierr);
177 ierr = PetscArrayzero(work,A->cmap->N);CHKERRQ(ierr);
178 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
179 for (i=rstart; i<rend; i++) {
180 ierr = MatGetRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr);
181 for (j=0; j<nnz; j++) {
182 work[cols[j]] += mvalues[j]*PetscConj(mvalues[j]);
183 }
184 ierr = MatRestoreRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr);
185 }
186 ierr = MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));CHKERRQ(ierr);
187 rstart = N->cmap->rstart;
188 rend = N->cmap->rend;
189 ierr = VecGetArray(v,&values);CHKERRQ(ierr);
190 ierr = PetscArraycpy(values,diag+rstart,rend-rstart);CHKERRQ(ierr);
191 ierr = VecRestoreArray(v,&values);CHKERRQ(ierr);
192 ierr = PetscFree2(diag,work);CHKERRQ(ierr);
193 ierr = VecScale(v,Na->scale);CHKERRQ(ierr);
194 PetscFunctionReturn(0);
195 }
196
197 /*@
198 MatCreateNormalHermitian - Creates a new matrix object that behaves like (A*)'*A.
199
200 Collective on Mat
201
202 Input Parameter:
203 . A - the (possibly rectangular complex) matrix
204
205 Output Parameter:
206 . N - the matrix that represents (A*)'*A
207
208 Level: intermediate
209
210 Notes:
211 The product (A*)'*A is NOT actually formed! Rather the new matrix
212 object performs the matrix-vector product by first multiplying by
213 A and then (A*)'
214 @*/
MatCreateNormalHermitian(Mat A,Mat * N)215 PetscErrorCode MatCreateNormalHermitian(Mat A,Mat *N)
216 {
217 PetscErrorCode ierr;
218 PetscInt m,n;
219 Mat_Normal *Na;
220
221 PetscFunctionBegin;
222 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
223 ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr);
224 ierr = MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
225 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATNORMALHERMITIAN);CHKERRQ(ierr);
226
227 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr);
228 (*N)->data = (void*) Na;
229 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
230 Na->A = A;
231 Na->scale = 1.0;
232
233 ierr = VecCreateMPI(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,&Na->w);CHKERRQ(ierr);
234
235 (*N)->ops->destroy = MatDestroyHermitian_Normal;
236 (*N)->ops->mult = MatMultHermitian_Normal;
237 (*N)->ops->multtranspose = MatMultHermitianTranspose_Normal;
238 (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_Normal;
239 (*N)->ops->multadd = MatMultHermitianAdd_Normal;
240 (*N)->ops->getdiagonal = MatGetDiagonalHermitian_Normal;
241 (*N)->ops->scale = MatScaleHermitian_Normal;
242 (*N)->ops->diagonalscale = MatDiagonalScaleHermitian_Normal;
243 (*N)->assembled = PETSC_TRUE;
244 (*N)->cmap->N = A->cmap->N;
245 (*N)->rmap->N = A->cmap->N;
246 (*N)->cmap->n = A->cmap->n;
247 (*N)->rmap->n = A->cmap->n;
248 PetscFunctionReturn(0);
249 }
250
251