1 #define PETSCMAT_DLL 2 3 #include "private/matimpl.h" /*I "petscmat.h" I*/ 4 5 typedef struct { 6 Mat A,left,right; /* left and right scaling not yet implemented */ 7 Vec w; 8 PetscScalar scale; 9 } Mat_Normal; 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "MatMult_Normal" 13 PetscErrorCode MatMult_Normal(Mat N,Vec x,Vec y) 14 { 15 Mat_Normal *Na = (Mat_Normal*)N->data; 16 PetscErrorCode ierr; 17 18 PetscFunctionBegin; 19 ierr = MatMult(Na->A,x,Na->w);CHKERRQ(ierr); 20 ierr = MatMultTranspose(Na->A,Na->w,y);CHKERRQ(ierr); 21 ierr = VecScale(y,Na->scale);CHKERRQ(ierr); 22 PetscFunctionReturn(0); 23 } 24 25 #undef __FUNCT__ 26 #define __FUNCT__ "MatMultAdd_Normal" 27 PetscErrorCode MatMultAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3) 28 { 29 Mat_Normal *Na = (Mat_Normal*)N->data; 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 ierr = MatMult(Na->A,v1,Na->w);CHKERRQ(ierr); 34 ierr = VecScale(Na->w,Na->scale);CHKERRQ(ierr); 35 ierr = MatMultTransposeAdd(Na->A,Na->w,v2,v3);CHKERRQ(ierr); 36 PetscFunctionReturn(0); 37 } 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "MatMultTranspose_Normal" 41 PetscErrorCode MatMultTranspose_Normal(Mat N,Vec x,Vec y) 42 { 43 Mat_Normal *Na = (Mat_Normal*)N->data; 44 PetscErrorCode ierr; 45 46 PetscFunctionBegin; 47 ierr = MatMult(Na->A,x,Na->w);CHKERRQ(ierr); 48 ierr = MatMultTranspose(Na->A,Na->w,y);CHKERRQ(ierr); 49 ierr = VecScale(y,Na->scale);CHKERRQ(ierr); 50 PetscFunctionReturn(0); 51 } 52 53 #undef __FUNCT__ 54 #define __FUNCT__ "MatMultTransposeAdd_Normal" 55 PetscErrorCode MatMultTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3) 56 { 57 Mat_Normal *Na = (Mat_Normal*)N->data; 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 ierr = MatMult(Na->A,v1,Na->w);CHKERRQ(ierr); 62 ierr = VecScale(Na->w,Na->scale);CHKERRQ(ierr); 63 ierr = MatMultTransposeAdd(Na->A,Na->w,v2,v3);CHKERRQ(ierr); 64 PetscFunctionReturn(0); 65 } 66 67 #undef __FUNCT__ 68 #define __FUNCT__ "MatDestroy_Normal" 69 PetscErrorCode MatDestroy_Normal(Mat N) 70 { 71 Mat_Normal *Na = (Mat_Normal*)N->data; 72 PetscErrorCode ierr; 73 74 PetscFunctionBegin; 75 if (Na->A) { ierr = MatDestroy(Na->A);CHKERRQ(ierr); } 76 if (Na->w) { ierr = VecDestroy(Na->w);CHKERRQ(ierr); } 77 ierr = PetscFree(Na);CHKERRQ(ierr); 78 PetscFunctionReturn(0); 79 } 80 81 /* 82 Slow, nonscalable version 83 */ 84 #undef __FUNCT__ 85 #define __FUNCT__ "MatGetDiagonal_Normal" 86 PetscErrorCode MatGetDiagonal_Normal(Mat N,Vec v) 87 { 88 Mat_Normal *Na = (Mat_Normal*)N->data; 89 Mat A = Na->A; 90 PetscErrorCode ierr; 91 PetscInt i,j,rstart,rend,nnz; 92 const PetscInt *cols; 93 PetscScalar *diag,*work,*values; 94 const PetscScalar *mvalues; 95 96 PetscFunctionBegin; 97 ierr = PetscMalloc(2*A->cmap->N*sizeof(PetscScalar),&diag);CHKERRQ(ierr); 98 work = diag + A->cmap->N; 99 ierr = PetscMemzero(work,A->cmap->N*sizeof(PetscScalar));CHKERRQ(ierr); 100 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 101 for (i=rstart; i<rend; i++) { 102 ierr = MatGetRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr); 103 for (j=0; j<nnz; j++) { 104 work[cols[j]] += mvalues[j]*mvalues[j]; 105 } 106 ierr = MatRestoreRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr); 107 } 108 ierr = MPI_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPI_SUM,((PetscObject)N)->comm);CHKERRQ(ierr); 109 rstart = N->cmap->rstart; 110 rend = N->cmap->rend; 111 ierr = VecGetArray(v,&values);CHKERRQ(ierr); 112 ierr = PetscMemcpy(values,diag+rstart,(rend-rstart)*sizeof(PetscScalar));CHKERRQ(ierr); 113 ierr = VecRestoreArray(v,&values);CHKERRQ(ierr); 114 ierr = PetscFree(diag);CHKERRQ(ierr); 115 ierr = VecScale(v,Na->scale);CHKERRQ(ierr); 116 PetscFunctionReturn(0); 117 } 118 119 #undef __FUNCT__ 120 #define __FUNCT__ "MatCreateNormal" 121 /*@ 122 MatCreateNormal - Creates a new matrix object that behaves like A'*A. 123 124 Collective on Mat 125 126 Input Parameter: 127 . A - the (possibly rectangular) matrix 128 129 Output Parameter: 130 . N - the matrix that represents A'*A 131 132 Level: intermediate 133 134 Notes: The product A'*A is NOT actually formed! Rather the new matrix 135 object performs the matrix-vector product by first multiplying by 136 A and then A' 137 @*/ 138 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateNormal(Mat A,Mat *N) 139 { 140 PetscErrorCode ierr; 141 PetscInt m,n; 142 Mat_Normal *Na; 143 144 PetscFunctionBegin; 145 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 146 ierr = MatCreate(((PetscObject)A)->comm,N);CHKERRQ(ierr); 147 ierr = MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 148 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATNORMAL);CHKERRQ(ierr); 149 150 ierr = PetscNewLog(*N,Mat_Normal,&Na);CHKERRQ(ierr); 151 (*N)->data = (void*) Na; 152 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 153 Na->A = A; 154 Na->scale = 1.0; 155 156 ierr = VecCreateMPI(((PetscObject)A)->comm,m,PETSC_DECIDE,&Na->w);CHKERRQ(ierr); 157 (*N)->ops->destroy = MatDestroy_Normal; 158 (*N)->ops->mult = MatMult_Normal; 159 (*N)->ops->multtranspose = MatMultTranspose_Normal; 160 (*N)->ops->multtransposeadd = MatMultTransposeAdd_Normal; 161 (*N)->ops->multadd = MatMultAdd_Normal; 162 (*N)->ops->getdiagonal = MatGetDiagonal_Normal; 163 (*N)->assembled = PETSC_TRUE; 164 (*N)->cmap->N = A->cmap->N; 165 (*N)->rmap->N = A->cmap->N; 166 (*N)->cmap->n = A->cmap->n; 167 (*N)->rmap->n = A->cmap->n; 168 PetscFunctionReturn(0); 169 } 170 171