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