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 PetscFunctionReturn(0); 160 } 161 162 /* 163 Slow, nonscalable version 164 */ 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 @*/ 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