1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 typedef struct { 5 Mat A; /* sparse matrix */ 6 Mat U,V; /* dense tall-skinny matrices */ 7 Vec c; /* sequential vector containing the diagonal of C */ 8 Vec work1,work2; /* sequential vectors that hold partial products */ 9 PetscMPIInt nwork; /* length of work vectors */ 10 } Mat_LRC; 11 12 13 PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat,Vec,Vec); 14 PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat,Vec,Vec,Vec); 15 PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat,Vec,Vec); 16 17 #undef __FUNCT__ 18 #define __FUNCT__ "MatMult_LRC" 19 PetscErrorCode MatMult_LRC(Mat N,Vec x,Vec y) 20 { 21 Mat_LRC *Na = (Mat_LRC*)N->data; 22 PetscErrorCode ierr; 23 PetscScalar *w1,*w2; 24 25 PetscFunctionBegin; 26 /* multiply the local part of V with the local part of x */ 27 /* note in this call x is treated as a sequential vector */ 28 ierr = MatMultTranspose_SeqDense(Na->V,x,Na->work1);CHKERRQ(ierr); 29 30 /* form the sum of all the local multiplies: this is work2 = V'*x = 31 sum_{all processors} work1 */ 32 ierr = VecGetArray(Na->work1,&w1);CHKERRQ(ierr); 33 ierr = VecGetArray(Na->work2,&w2);CHKERRQ(ierr); 34 ierr = MPIU_Allreduce(w1,w2,Na->nwork,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));CHKERRQ(ierr); 35 ierr = VecRestoreArray(Na->work1,&w1);CHKERRQ(ierr); 36 ierr = VecRestoreArray(Na->work2,&w2);CHKERRQ(ierr); 37 38 if (Na->c) { /* work2 = C*work2 */ 39 ierr = VecPointwiseMult(Na->work2,Na->c,Na->work2);CHKERRQ(ierr); 40 } 41 42 if (Na->A) { 43 /* form y = A*x */ 44 ierr = MatMult(Na->A,x,y);CHKERRQ(ierr); 45 /* multiply-add y = y + U*work2 */ 46 /* note in this call y is treated as a sequential vector */ 47 ierr = MatMultAdd_SeqDense(Na->U,Na->work2,y,y);CHKERRQ(ierr); 48 } else { 49 /* multiply y = U*work2 */ 50 ierr = MatMult_SeqDense(Na->U,Na->work2,y);CHKERRQ(ierr); 51 } 52 PetscFunctionReturn(0); 53 } 54 55 #undef __FUNCT__ 56 #define __FUNCT__ "MatDestroy_LRC" 57 PetscErrorCode MatDestroy_LRC(Mat N) 58 { 59 Mat_LRC *Na = (Mat_LRC*)N->data; 60 PetscErrorCode ierr; 61 62 PetscFunctionBegin; 63 ierr = MatDestroy(&Na->A);CHKERRQ(ierr); 64 ierr = MatDestroy(&Na->U);CHKERRQ(ierr); 65 ierr = MatDestroy(&Na->V);CHKERRQ(ierr); 66 ierr = VecDestroy(&Na->c);CHKERRQ(ierr); 67 ierr = VecDestroy(&Na->work1);CHKERRQ(ierr); 68 ierr = VecDestroy(&Na->work2);CHKERRQ(ierr); 69 ierr = PetscFree(N->data);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "MatCreateLRC" 75 /*@ 76 MatCreateLRC - Creates a new matrix object that behaves like A + U*C*V' 77 78 Collective on Mat 79 80 Input Parameters: 81 + A - the (sparse) matrix (can be NULL) 82 . U, V - two dense rectangular (tall and skinny) matrices 83 - c - a sequential vector containing the diagonal of C (can be NULL) 84 85 Output Parameter: 86 . N - the matrix that represents A + U*C*V' 87 88 Notes: 89 The matrix A + U*C*V' is not formed! Rather the new matrix 90 object performs the matrix-vector product by first multiplying by 91 A and then adding the other term. 92 93 C is a diagonal matrix (represented as a vector) of order k, 94 where k is the number of columns of both U and V. 95 96 If A is NULL then the new object behaves like a low-rank matrix U*C*V'. 97 98 Use V=U (or V=NULL) for a symmetric low-rank correction, A + U*C*U'. 99 100 If c is NULL then the low-rank correction is just U*V'. 101 102 Level: intermediate 103 @*/ 104 PetscErrorCode MatCreateLRC(Mat A,Mat U,Vec c,Mat V,Mat *N) 105 { 106 PetscErrorCode ierr; 107 PetscBool match; 108 PetscInt m,n,k,m1,n1,k1; 109 Mat_LRC *Na; 110 111 PetscFunctionBegin; 112 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,1); 113 PetscValidHeaderSpecific(U,MAT_CLASSID,2); 114 if (c) PetscValidHeaderSpecific(c,VEC_CLASSID,3); 115 if (V) PetscValidHeaderSpecific(V,MAT_CLASSID,4); 116 else V=U; 117 if (A) PetscCheckSameComm(A,1,U,2); 118 PetscCheckSameComm(U,2,V,4); 119 120 ierr = PetscObjectTypeCompareAny((PetscObject)U,&match,MATSEQDENSE,MATMPIDENSE,"");CHKERRQ(ierr); 121 if (!match) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_SUP,"Matrix U must be of type dense"); 122 if (V) { 123 ierr = PetscObjectTypeCompareAny((PetscObject)V,&match,MATSEQDENSE,MATMPIDENSE,"");CHKERRQ(ierr); 124 if (!match) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Matrix V must be of type dense"); 125 } 126 127 ierr = MatGetSize(U,NULL,&k);CHKERRQ(ierr); 128 ierr = MatGetSize(V,NULL,&k1);CHKERRQ(ierr); 129 if (k!=k1) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_INCOMP,"U and V have different number of columns"); 130 ierr = MatGetLocalSize(U,&m,NULL);CHKERRQ(ierr); 131 ierr = MatGetLocalSize(V,&n,NULL);CHKERRQ(ierr); 132 if (A) { 133 ierr = MatGetLocalSize(A,&m1,&n1);CHKERRQ(ierr); 134 if (m!=m1) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_INCOMP,"Local dimensions of U and A do not match"); 135 if (n!=n1) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Local dimensions of V and A do not match"); 136 } 137 if (c) { 138 ierr = VecGetSize(c,&k1);CHKERRQ(ierr); 139 if (k!=k1) SETERRQ(PetscObjectComm((PetscObject)c),PETSC_ERR_ARG_INCOMP,"The length of c does not match the number of columns of U and V"); 140 ierr = VecGetLocalSize(c,&k1);CHKERRQ(ierr); 141 if (k!=k1) SETERRQ(PetscObjectComm((PetscObject)c),PETSC_ERR_ARG_INCOMP,"c must be a sequential vector"); 142 } 143 144 ierr = MatCreate(PetscObjectComm((PetscObject)U),N);CHKERRQ(ierr); 145 ierr = MatSetSizes(*N,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 146 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATLRC);CHKERRQ(ierr); 147 148 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr); 149 (*N)->data = (void*)Na; 150 Na->A = A; 151 Na->c = c; 152 153 ierr = MatDenseGetLocalMatrix(U,&Na->U);CHKERRQ(ierr); 154 ierr = MatDenseGetLocalMatrix(V,&Na->V);CHKERRQ(ierr); 155 if (A) { ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); } 156 ierr = PetscObjectReference((PetscObject)Na->U);CHKERRQ(ierr); 157 ierr = PetscObjectReference((PetscObject)Na->V);CHKERRQ(ierr); 158 if (c) { ierr = PetscObjectReference((PetscObject)c);CHKERRQ(ierr); } 159 160 ierr = VecCreateSeq(PETSC_COMM_SELF,U->cmap->N,&Na->work1);CHKERRQ(ierr); 161 ierr = VecDuplicate(Na->work1,&Na->work2);CHKERRQ(ierr); 162 ierr = PetscMPIIntCast(U->cmap->N,&Na->nwork);CHKERRQ(ierr); 163 164 (*N)->ops->destroy = MatDestroy_LRC; 165 (*N)->ops->mult = MatMult_LRC; 166 (*N)->assembled = PETSC_TRUE; 167 (*N)->preallocated = PETSC_TRUE; 168 (*N)->cmap->N = V->rmap->N; 169 (*N)->rmap->N = U->rmap->N; 170 (*N)->cmap->n = V->rmap->n; 171 (*N)->rmap->n = U->rmap->n; 172 PetscFunctionReturn(0); 173 } 174 175