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 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 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 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 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 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 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 ierr = PetscObjectComposeFunction((PetscObject)N,"MatNormalGetMatHermitian_C",NULL);CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 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; 175 176 PetscFunctionBegin; 177 ierr = PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);CHKERRQ(ierr); 178 ierr = PetscArrayzero(work,A->cmap->N);CHKERRQ(ierr); 179 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 180 for (i=rstart; i<rend; i++) { 181 ierr = MatGetRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr); 182 for (j=0; j<nnz; j++) { 183 work[cols[j]] += mvalues[j]*PetscConj(mvalues[j]); 184 } 185 ierr = MatRestoreRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr); 186 } 187 ierr = MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));CHKERRMPI(ierr); 188 rstart = N->cmap->rstart; 189 rend = N->cmap->rend; 190 ierr = VecGetArray(v,&values);CHKERRQ(ierr); 191 ierr = PetscArraycpy(values,diag+rstart,rend-rstart);CHKERRQ(ierr); 192 ierr = VecRestoreArray(v,&values);CHKERRQ(ierr); 193 ierr = PetscFree2(diag,work);CHKERRQ(ierr); 194 ierr = VecScale(v,Na->scale);CHKERRQ(ierr); 195 PetscFunctionReturn(0); 196 } 197 198 PetscErrorCode MatNormalGetMatHermitian_Normal(Mat A,Mat *M) 199 { 200 Mat_Normal *Aa = (Mat_Normal*)A->data; 201 202 PetscFunctionBegin; 203 *M = Aa->A; 204 PetscFunctionReturn(0); 205 } 206 207 /*@ 208 MatNormalHermitianGetMat - Gets the Mat object stored inside a MATNORMALHERMITIAN 209 210 Logically collective on Mat 211 212 Input Parameter: 213 . A - the MATNORMALHERMITIAN matrix 214 215 Output Parameter: 216 . M - the matrix object stored inside A 217 218 Level: intermediate 219 220 .seealso: MatCreateNormalHermitian() 221 222 @*/ 223 PetscErrorCode MatNormalHermitianGetMat(Mat A,Mat *M) 224 { 225 PetscErrorCode ierr; 226 227 PetscFunctionBegin; 228 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 229 PetscValidType(A,1); 230 PetscValidPointer(M,2); 231 ierr = PetscUseMethod(A,"MatNormalGetMatHermitian_C",(Mat,Mat*),(A,M));CHKERRQ(ierr); 232 PetscFunctionReturn(0); 233 } 234 235 /*@ 236 MatCreateNormalHermitian - Creates a new matrix object that behaves like (A*)'*A. 237 238 Collective on Mat 239 240 Input Parameter: 241 . A - the (possibly rectangular complex) matrix 242 243 Output Parameter: 244 . N - the matrix that represents (A*)'*A 245 246 Level: intermediate 247 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 { 255 PetscErrorCode ierr; 256 PetscInt m,n; 257 Mat_Normal *Na; 258 VecType vtype; 259 260 PetscFunctionBegin; 261 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 262 ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr); 263 ierr = MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 264 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATNORMALHERMITIAN);CHKERRQ(ierr); 265 ierr = PetscLayoutReference(A->cmap,&(*N)->rmap);CHKERRQ(ierr); 266 ierr = PetscLayoutReference(A->cmap,&(*N)->cmap);CHKERRQ(ierr); 267 268 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr); 269 (*N)->data = (void*) Na; 270 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 271 Na->A = A; 272 Na->scale = 1.0; 273 274 ierr = MatCreateVecs(A,NULL,&Na->w);CHKERRQ(ierr); 275 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; 286 287 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatNormalGetMatHermitian_C",MatNormalGetMatHermitian_Normal);CHKERRQ(ierr); 288 ierr = MatSetOption(*N,MAT_HERMITIAN,PETSC_TRUE);CHKERRQ(ierr); 289 ierr = MatGetVecType(A,&vtype);CHKERRQ(ierr); 290 ierr = MatSetVecType(*N,vtype);CHKERRQ(ierr); 291 #if defined(PETSC_HAVE_DEVICE) 292 ierr = MatBindToCPU(*N,A->boundtocpu);CHKERRQ(ierr); 293 #endif 294 PetscFunctionReturn(0); 295 } 296 297