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