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