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