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 vectors that hold partial products */ 7 PetscMPIInt nwork; /* length of work vectors */ 8 } Mat_LRC; 9 10 11 PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat,Vec,Vec); 12 PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat,Vec,Vec,Vec); 13 PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat,Vec,Vec); 14 15 #undef __FUNCT__ 16 #define __FUNCT__ "MatMult_LRC" 17 PetscErrorCode MatMult_LRC(Mat N,Vec x,Vec y) 18 { 19 Mat_LRC *Na = (Mat_LRC*)N->data; 20 PetscErrorCode ierr; 21 PetscScalar *w1,*w2; 22 23 PetscFunctionBegin; 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 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 if (Na->A) { 37 /* form y = A*x */ 38 ierr = MatMult(Na->A,x,y);CHKERRQ(ierr); 39 /* multiply-add y = y + U*work2 */ 40 /* note in this call y is treated as a sequential vector */ 41 ierr = MatMultAdd_SeqDense(Na->U,Na->work2,y,y);CHKERRQ(ierr); 42 } else { 43 /* multiply y = U*work2 */ 44 ierr = MatMult_SeqDense(Na->U,Na->work2,y);CHKERRQ(ierr); 45 } 46 PetscFunctionReturn(0); 47 } 48 49 #undef __FUNCT__ 50 #define __FUNCT__ "MatDestroy_LRC" 51 PetscErrorCode MatDestroy_LRC(Mat N) 52 { 53 Mat_LRC *Na = (Mat_LRC*)N->data; 54 PetscErrorCode ierr; 55 56 PetscFunctionBegin; 57 ierr = MatDestroy(&Na->A);CHKERRQ(ierr); 58 ierr = MatDestroy(&Na->U);CHKERRQ(ierr); 59 ierr = MatDestroy(&Na->V);CHKERRQ(ierr); 60 ierr = VecDestroy(&Na->work1);CHKERRQ(ierr); 61 ierr = VecDestroy(&Na->work2);CHKERRQ(ierr); 62 ierr = PetscFree(N->data);CHKERRQ(ierr); 63 PetscFunctionReturn(0); 64 } 65 66 #undef __FUNCT__ 67 #define __FUNCT__ "MatCreateLRC" 68 /*@ 69 MatCreateLRC - Creates a new matrix object that behaves like A + U*V' 70 71 Collective on Mat 72 73 Input Parameters: 74 + A - the (sparse) matrix (can be NULL) 75 - U, V - two dense rectangular (tall and skinny) matrices 76 77 Output Parameter: 78 . N - the matrix that represents A + U*V' 79 80 Notes: 81 The matrix A + U*V' is not formed! Rather the new matrix 82 object performs the matrix-vector product by first multiplying by 83 A and then adding the other term. 84 85 If A is NULL then the new object behaves like a low-rank matrix U*V'. 86 87 Use V=U (or V=NULL) for a symmetric low-rank correction, A + U*U'. 88 89 Level: intermediate 90 @*/ 91 PetscErrorCode MatCreateLRC(Mat A,Mat U,Mat V,Mat *N) 92 { 93 PetscErrorCode ierr; 94 PetscInt m,n,k,m1,n1,k1; 95 Mat_LRC *Na; 96 97 PetscFunctionBegin; 98 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,1); 99 PetscValidHeaderSpecific(U,MAT_CLASSID,2); 100 if (V) PetscValidHeaderSpecific(V,MAT_CLASSID,3); 101 else V=U; 102 if (A) PetscCheckSameComm(A,1,U,2); 103 PetscCheckSameComm(U,2,V,3); 104 105 ierr = MatGetSize(U,NULL,&k);CHKERRQ(ierr); 106 ierr = MatGetSize(V,NULL,&k1);CHKERRQ(ierr); 107 if (k!=k1) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_INCOMP,"U and V have different number of columns"); 108 ierr = MatGetLocalSize(U,&m,NULL);CHKERRQ(ierr); 109 ierr = MatGetLocalSize(V,&n,NULL);CHKERRQ(ierr); 110 if (A) { 111 ierr = MatGetLocalSize(A,&m1,&n1);CHKERRQ(ierr); 112 if (m!=m1) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_INCOMP,"Local dimensions of U and A do not match"); 113 if (n!=n1) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Local dimensions of V and A do not match"); 114 } 115 116 ierr = MatCreate(PetscObjectComm((PetscObject)U),N);CHKERRQ(ierr); 117 ierr = MatSetSizes(*N,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 118 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATLRC);CHKERRQ(ierr); 119 120 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr); 121 (*N)->data = (void*)Na; 122 Na->A = A; 123 124 ierr = MatDenseGetLocalMatrix(U,&Na->U);CHKERRQ(ierr); 125 ierr = MatDenseGetLocalMatrix(V,&Na->V);CHKERRQ(ierr); 126 if (A) { ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); } 127 ierr = PetscObjectReference((PetscObject)Na->U);CHKERRQ(ierr); 128 ierr = PetscObjectReference((PetscObject)Na->V);CHKERRQ(ierr); 129 130 ierr = VecCreateSeq(PETSC_COMM_SELF,U->cmap->N,&Na->work1);CHKERRQ(ierr); 131 ierr = VecDuplicate(Na->work1,&Na->work2);CHKERRQ(ierr); 132 Na->nwork = U->cmap->N; 133 134 (*N)->ops->destroy = MatDestroy_LRC; 135 (*N)->ops->mult = MatMult_LRC; 136 (*N)->assembled = PETSC_TRUE; 137 (*N)->cmap->N = V->rmap->N; 138 (*N)->rmap->N = U->rmap->N; 139 (*N)->cmap->n = V->rmap->n; 140 (*N)->rmap->n = U->rmap->n; 141 PetscFunctionReturn(0); 142 } 143 144