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