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