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