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