1 #define PETSCMAT_DLL 2 3 #include "src/mat/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__ "MatDestroy_Normal" 25 PetscErrorCode MatDestroy_Normal(Mat N) 26 { 27 Mat_Normal *Na = (Mat_Normal*)N->data; 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 ierr = PetscObjectDereference((PetscObject)Na->A);CHKERRQ(ierr); 32 ierr = VecDestroy(Na->w);CHKERRQ(ierr); 33 ierr = PetscFree(Na);CHKERRQ(ierr); 34 PetscFunctionReturn(0); 35 } 36 37 /* 38 Slow, nonscalable version 39 */ 40 #undef __FUNCT__ 41 #define __FUNCT__ "MatGetDiagonal_Normal" 42 PetscErrorCode MatGetDiagonal_Normal(Mat N,Vec v) 43 { 44 Mat_Normal *Na = (Mat_Normal*)N->data; 45 Mat A = Na->A; 46 PetscErrorCode ierr; 47 PetscInt i,j,rstart,rend,nnz; 48 const PetscInt *cols; 49 PetscScalar *diag,*work,*values; 50 const PetscScalar *mvalues; 51 PetscMap cmap; 52 53 PetscFunctionBegin; 54 ierr = PetscMalloc(2*A->N*sizeof(PetscScalar),&diag);CHKERRQ(ierr); 55 work = diag + A->N; 56 ierr = PetscMemzero(work,A->N*sizeof(PetscScalar));CHKERRQ(ierr); 57 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 58 for (i=rstart; i<rend; i++) { 59 ierr = MatGetRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr); 60 for (j=0; j<nnz; j++) { 61 work[cols[j]] += mvalues[j]*mvalues[j]; 62 } 63 ierr = MatRestoreRow(A,i,&nnz,&cols,&mvalues);CHKERRQ(ierr); 64 } 65 ierr = MPI_Allreduce(work,diag,A->N,MPIU_SCALAR,MPI_SUM,N->comm);CHKERRQ(ierr); 66 ierr = MatGetPetscMaps(A,PETSC_NULL,&cmap);CHKERRQ(ierr); 67 ierr = PetscMapGetLocalRange(cmap,&rstart,&rend);CHKERRQ(ierr); 68 ierr = VecGetArray(v,&values);CHKERRQ(ierr); 69 ierr = PetscMemcpy(values,diag+rstart,(rend-rstart)*sizeof(PetscScalar));CHKERRQ(ierr); 70 ierr = VecRestoreArray(v,&values);CHKERRQ(ierr); 71 ierr = PetscFree(diag);CHKERRQ(ierr); 72 PetscFunctionReturn(0); 73 } 74 75 #undef __FUNCT__ 76 #define __FUNCT__ "MatCreateNormal" 77 /*@ 78 MatCreateNormal - Creates a new matrix object that behaves like A'*A. 79 80 Collective on Mat 81 82 Input Parameter: 83 . A - the (possibly rectangular) matrix 84 85 Output Parameter: 86 . N - the matrix that represents A'*A 87 88 Level: intermediate 89 90 Notes: The product A'*A is NOT actually formed! Rather the new matrix 91 object performs the matrix-vector product by first multiplying by 92 A and then A' 93 @*/ 94 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateNormal(Mat A,Mat *N) 95 { 96 PetscErrorCode ierr; 97 PetscInt m,n; 98 Mat_Normal *Na; 99 100 PetscFunctionBegin; 101 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 102 ierr = MatCreate(A->comm,n,n,PETSC_DECIDE,PETSC_DECIDE,N);CHKERRQ(ierr); 103 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATNORMAL);CHKERRQ(ierr); 104 105 ierr = PetscNew(Mat_Normal,&Na);CHKERRQ(ierr); 106 Na->A = A; 107 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 108 (*N)->data = (void*) Na; 109 110 ierr = VecCreateMPI(A->comm,m,PETSC_DECIDE,&Na->w);CHKERRQ(ierr); 111 (*N)->ops->destroy = MatDestroy_Normal; 112 (*N)->ops->mult = MatMult_Normal; 113 (*N)->ops->getdiagonal = MatGetDiagonal_Normal; 114 (*N)->assembled = PETSC_TRUE; 115 (*N)->N = A->N; 116 (*N)->M = A->N; 117 (*N)->n = A->n; 118 (*N)->m = A->n; 119 PetscFunctionReturn(0); 120 } 121 122