Actual source code: lrc.c
2: #include <petsc/private/matimpl.h>
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: Vec xl,yl; /* auxiliary sequential vectors for matmult operation */
11: } Mat_LRC;
13: PetscErrorCode MatMult_LRC(Mat N,Vec x,Vec y)
14: {
15: Mat_LRC *Na = (Mat_LRC*)N->data;
16: Mat Uloc,Vloc;
17: PetscErrorCode ierr;
18: PetscScalar *w1,*w2;
19: const PetscScalar *a;
22: VecGetArrayRead(x,&a);
23: VecPlaceArray(Na->xl,a);
24: VecGetLocalVector(y,Na->yl);
25: MatDenseGetLocalMatrix(Na->U,&Uloc);
26: MatDenseGetLocalMatrix(Na->V,&Vloc);
28: /* multiply the local part of V with the local part of x */
29: #if defined(PETSC_USE_COMPLEX)
30: MatMultHermitianTranspose(Vloc,Na->xl,Na->work1);
31: #else
32: MatMultTranspose(Vloc,Na->xl,Na->work1);
33: #endif
35: /* form the sum of all the local multiplies: this is work2 = V'*x =
36: sum_{all processors} work1 */
37: VecGetArray(Na->work1,&w1);
38: VecGetArray(Na->work2,&w2);
39: MPIU_Allreduce(w1,w2,Na->nwork,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));
40: VecRestoreArray(Na->work1,&w1);
41: VecRestoreArray(Na->work2,&w2);
43: if (Na->c) { /* work2 = C*work2 */
44: VecPointwiseMult(Na->work2,Na->c,Na->work2);
45: }
47: if (Na->A) {
48: /* form y = A*x */
49: MatMult(Na->A,x,y);
50: /* multiply-add y = y + U*work2 */
51: MatMultAdd(Uloc,Na->work2,Na->yl,Na->yl);
52: } else {
53: /* multiply y = U*work2 */
54: MatMult(Uloc,Na->work2,Na->yl);
55: }
57: VecRestoreArrayRead(x,&a);
58: VecResetArray(Na->xl);
59: VecRestoreLocalVector(y,Na->yl);
60: return(0);
61: }
63: PetscErrorCode MatDestroy_LRC(Mat N)
64: {
65: Mat_LRC *Na = (Mat_LRC*)N->data;
69: MatDestroy(&Na->A);
70: MatDestroy(&Na->U);
71: MatDestroy(&Na->V);
72: VecDestroy(&Na->c);
73: VecDestroy(&Na->work1);
74: VecDestroy(&Na->work2);
75: VecDestroy(&Na->xl);
76: VecDestroy(&Na->yl);
77: PetscFree(N->data);
78: PetscObjectComposeFunction((PetscObject)N,"MatLRCGetMats_C",0);
79: return(0);
80: }
82: PetscErrorCode MatLRCGetMats_LRC(Mat N,Mat *A,Mat *U,Vec *c,Mat *V)
83: {
84: Mat_LRC *Na = (Mat_LRC*)N->data;
87: if (A) *A = Na->A;
88: if (U) *U = Na->U;
89: if (c) *c = Na->c;
90: if (V) *V = Na->V;
91: return(0);
92: }
94: /*@
95: MatLRCGetMats - Returns the constituents of an LRC matrix
97: Collective on Mat
99: Input Parameter:
100: . N - matrix of type LRC
102: Output Parameters:
103: + A - the (sparse) matrix
104: . U - first dense rectangular (tall and skinny) matrix
105: . c - a sequential vector containing the diagonal of C
106: - V - second dense rectangular (tall and skinny) matrix
108: Note:
109: The returned matrices need not be destroyed by the caller.
111: Level: intermediate
113: .seealso: MatCreateLRC()
114: @*/
115: PetscErrorCode MatLRCGetMats(Mat N,Mat *A,Mat *U,Vec *c,Mat *V)
116: {
120: PetscUseMethod(N,"MatLRCGetMats_C",(Mat,Mat*,Mat*,Vec*,Mat*),(N,A,U,c,V));
121: return(0);
122: }
124: /*@
125: MatCreateLRC - Creates a new matrix object that behaves like A + U*C*V'
127: Collective on Mat
129: Input Parameters:
130: + A - the (sparse) matrix (can be NULL)
131: . U, V - two dense rectangular (tall and skinny) matrices
132: - c - a sequential vector containing the diagonal of C (can be NULL)
134: Output Parameter:
135: . N - the matrix that represents A + U*C*V'
137: Notes:
138: The matrix A + U*C*V' is not formed! Rather the new matrix
139: object performs the matrix-vector product by first multiplying by
140: A and then adding the other term.
142: C is a diagonal matrix (represented as a vector) of order k,
143: where k is the number of columns of both U and V.
145: If A is NULL then the new object behaves like a low-rank matrix U*C*V'.
147: Use V=U (or V=NULL) for a symmetric low-rank correction, A + U*C*U'.
149: If c is NULL then the low-rank correction is just U*V'.
151: Level: intermediate
153: .seealso: MatLRCGetMats()
154: @*/
155: PetscErrorCode MatCreateLRC(Mat A,Mat U,Vec c,Mat V,Mat *N)
156: {
158: PetscBool match;
159: PetscInt m,n,k,m1,n1,k1;
160: Mat_LRC *Na;
167: else V=U;
171: PetscObjectTypeCompareAny((PetscObject)U,&match,MATSEQDENSE,MATMPIDENSE,"");
172: if (!match) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_SUP,"Matrix U must be of type dense");
173: if (V) {
174: PetscObjectTypeCompareAny((PetscObject)V,&match,MATSEQDENSE,MATMPIDENSE,"");
175: if (!match) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Matrix V must be of type dense");
176: }
178: MatGetSize(U,NULL,&k);
179: MatGetSize(V,NULL,&k1);
180: if (k!=k1) SETERRQ2(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_INCOMP,"U and V have different number of columns (%D vs %D)",k,k1);
181: MatGetLocalSize(U,&m,NULL);
182: MatGetLocalSize(V,&n,NULL);
183: if (A) {
184: MatGetLocalSize(A,&m1,&n1);
185: if (m!=m1) SETERRQ2(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_INCOMP,"Local dimensions of U %D and A %D do not match",m,m1);
186: if (n!=n1) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Local dimensions of V %D and A %D do not match",n,n1);
187: }
188: if (c) {
189: VecGetSize(c,&k1);
190: if (k!=k1) SETERRQ2(PetscObjectComm((PetscObject)c),PETSC_ERR_ARG_INCOMP,"The length of c %D does not match the number of columns of U and V (%D)",k1,k);
191: VecGetLocalSize(c,&k1);
192: if (k!=k1) SETERRQ(PetscObjectComm((PetscObject)c),PETSC_ERR_ARG_INCOMP,"c must be a sequential vector");
193: }
195: MatCreate(PetscObjectComm((PetscObject)U),N);
196: MatSetSizes(*N,m,n,PETSC_DECIDE,PETSC_DECIDE);
197: PetscObjectChangeTypeName((PetscObject)*N,MATLRC);
199: PetscNewLog(*N,&Na);
200: (*N)->data = (void*)Na;
201: Na->A = A;
202: Na->U = U;
203: Na->c = c;
204: Na->V = V;
206: if (A) { PetscObjectReference((PetscObject)A); }
207: PetscObjectReference((PetscObject)Na->U);
208: PetscObjectReference((PetscObject)Na->V);
209: if (c) { PetscObjectReference((PetscObject)c); }
211: VecCreateSeq(PETSC_COMM_SELF,U->cmap->N,&Na->work1);
212: VecDuplicate(Na->work1,&Na->work2);
213: PetscMPIIntCast(U->cmap->N,&Na->nwork);
215: VecCreateSeqWithArray(PETSC_COMM_SELF,1,V->rmap->n,NULL,&Na->xl);
216: VecCreateSeqWithArray(PETSC_COMM_SELF,1,U->rmap->n,NULL,&Na->yl);
218: (*N)->ops->destroy = MatDestroy_LRC;
219: (*N)->ops->mult = MatMult_LRC;
220: (*N)->assembled = PETSC_TRUE;
221: (*N)->preallocated = PETSC_TRUE;
222: (*N)->cmap->N = V->rmap->N;
223: (*N)->rmap->N = U->rmap->N;
224: (*N)->cmap->n = V->rmap->n;
225: (*N)->rmap->n = U->rmap->n;
227: PetscObjectComposeFunction((PetscObject)(*N),"MatLRCGetMats_C",MatLRCGetMats_LRC);
228: return(0);
229: }