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