1 #define PETSCMAT_DLL 2 3 #include "include/private/matimpl.h" /*I "petscmat.h" I*/ 4 #include "src/mat/impls/dense/seq/dense.h" 5 6 typedef struct { 7 Mat A,U,V; 8 Vec work1,work2;/* Sequential (big) vectors that hold partial products */ 9 PetscMPIInt nwork; /* length of work vectors */ 10 } Mat_LRC; 11 12 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "MatMult_LRC" 16 PetscErrorCode MatMult_LRC(Mat N,Vec x,Vec y) 17 { 18 Mat_LRC *Na = (Mat_LRC*)N->data; 19 PetscErrorCode ierr; 20 PetscScalar *w1,*w2; 21 22 PetscFunctionBegin; 23 ierr = MatMult(Na->A,x,y);CHKERRQ(ierr); 24 25 /* multiply the local part of V with the local part of x */ 26 /* note in this call x is treated as a sequential vector */ 27 ierr = MatMultTranspose_SeqDense(Na->V,x,Na->work1);CHKERRQ(ierr); 28 29 /* Form the sum of all the local multiplies : this is work2 = V'*x = 30 sum_{all processors} work1 */ 31 32 ierr = VecGetArray(Na->work1,&w1);CHKERRQ(ierr); 33 ierr = VecGetArray(Na->work2,&w2);CHKERRQ(ierr); 34 ierr = MPI_Allreduce(w1,w2,Na->nwork,MPIU_SCALAR,PetscSum_Op,N->comm);CHKERRQ(ierr); 35 ierr = VecRestoreArray(Na->work1,&w1);CHKERRQ(ierr); 36 ierr = VecRestoreArray(Na->work2,&w2);CHKERRQ(ierr); 37 38 /* multiply-sub y = y + U*work2 */ 39 /* note in this call y is treated as a sequential vector */ 40 ierr = MatMultAdd_SeqDense(Na->U,Na->work2,y,y);CHKERRQ(ierr); 41 PetscFunctionReturn(0); 42 } 43 44 #undef __FUNCT__ 45 #define __FUNCT__ "MatDestroy_LRC" 46 PetscErrorCode MatDestroy_LRC(Mat N) 47 { 48 Mat_LRC *Na = (Mat_LRC*)N->data; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 ierr = PetscObjectDereference((PetscObject)Na->A);CHKERRQ(ierr); 53 ierr = PetscObjectDereference((PetscObject)Na->U);CHKERRQ(ierr); 54 ierr = PetscObjectDereference((PetscObject)Na->V);CHKERRQ(ierr); 55 ierr = VecDestroy(Na->work1);CHKERRQ(ierr); 56 ierr = VecDestroy(Na->work2);CHKERRQ(ierr); 57 ierr = PetscFree(Na);CHKERRQ(ierr); 58 PetscFunctionReturn(0); 59 } 60 61 62 #undef __FUNCT__ 63 #define __FUNCT__ "MatCreateLRC" 64 /*@ 65 MatCreateLRC - Creates a new matrix object that behaves like A + U*V' 66 67 Collective on Mat 68 69 Input Parameter: 70 + A - the (sparse) matrix 71 - U. V - two dense rectangular (tall and skinny) matrices 72 73 Output Parameter: 74 . N - the matrix that represents A + U*V' 75 76 Level: intermediate 77 78 Notes: The matrix A + U*V' formed! Rather the new matrix 79 object performs the matrix-vector product by first multiplying by 80 A and then adding the other term 81 @*/ 82 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateLRC(Mat A,Mat U, Mat V,Mat *N) 83 { 84 PetscErrorCode ierr; 85 PetscInt m,n; 86 Mat_LRC *Na; 87 88 PetscFunctionBegin; 89 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 90 ierr = MatCreate(A->comm,N);CHKERRQ(ierr); 91 ierr = MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 92 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATLRC);CHKERRQ(ierr); 93 94 ierr = PetscNewLog(*N,Mat_LRC,&Na);CHKERRQ(ierr); 95 (*N)->data = (void*) Na; 96 Na->A = A; 97 98 ierr = MatDenseGetLocalMatrix(U,&Na->U);CHKERRQ(ierr); 99 ierr = MatDenseGetLocalMatrix(V,&Na->V);CHKERRQ(ierr); 100 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 101 ierr = PetscObjectReference((PetscObject)Na->U);CHKERRQ(ierr); 102 ierr = PetscObjectReference((PetscObject)Na->V);CHKERRQ(ierr); 103 104 ierr = VecCreateSeq(PETSC_COMM_SELF,U->cmap.N,&Na->work1);CHKERRQ(ierr); 105 ierr = VecDuplicate(Na->work1,&Na->work2);CHKERRQ(ierr); 106 Na->nwork = U->cmap.N; 107 108 (*N)->ops->destroy = MatDestroy_LRC; 109 (*N)->ops->mult = MatMult_LRC; 110 (*N)->assembled = PETSC_TRUE; 111 (*N)->cmap.N = A->cmap.N; 112 (*N)->rmap.N = A->cmap.N; 113 (*N)->cmap.n = A->cmap.n; 114 (*N)->rmap.n = A->cmap.n; 115 PetscFunctionReturn(0); 116 } 117 118