Actual source code: dsgsvd.c
slepc-3.19.2 2023-09-05
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: #include <slepc/private/dsimpl.h>
11: #include <slepcblaslapack.h>
13: typedef struct {
14: PetscInt m; /* number of columns */
15: PetscInt p; /* number of rows of B */
16: PetscInt tm; /* number of rows of X after truncating */
17: PetscInt tp; /* number of rows of V after truncating */
18: } DS_GSVD;
20: PetscErrorCode DSAllocate_GSVD(DS ds,PetscInt ld)
21: {
22: PetscFunctionBegin;
23: PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
24: PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
25: PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
26: PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
27: PetscCall(DSAllocateMat_Private(ds,DS_MAT_V));
28: PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
29: PetscCall(DSAllocateMat_Private(ds,DS_MAT_D));
30: PetscCall(PetscFree(ds->perm));
31: PetscCall(PetscMalloc1(ld,&ds->perm));
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: /*
36: In compact form, A is either in form (a) or (b):
38: (a) (b)
39: lower bidiagonal with upper arrow (n=m+1) square upper bidiagonal with upper arrow (n=m)
40: 0 l k m-1
41: ----------------------------------------- 0 l k m-1
42: |* . | -----------------------------------------
43: | * . | |* . |
44: | * . | | * . |
45: | * . | | * . |
46: l |. . . . o o | l |. . . o o |
47: | o o | | o o |
48: | o o | | o o |
49: | o o | | o o |
50: | o o | | o o |
51: | o o | | o o |
52: k |. . . . . . . . . . o | k |. . . . . . . . . o x |
53: | x x | | x x |
54: | x x | | x x |
55: | x x | | x x |
56: | x x | | x x |
57: | x x | | x x |
58: | x x | | x x |
59: | x x | | x x |
60: | x x | | x x |
61: | x x| | x x|
62: n-1 | x| n-1 | x|
63: ----------------------------------------- -----------------------------------------
65: and B is square bidiagonal with upper arrow (p=m)
67: 0 l k m-1
68: -----------------------------------------
69: |* . |
70: | * . |
71: | * . |
72: | * . |
73: l |. . . . o o |
74: | o o |
75: | o o |
76: | o o |
77: | o o |
78: | o o |
79: k |. . . . . . . . . . o x |
80: | x x |
81: | x x |
82: | x x |
83: | x x |
84: | x x |
85: | x x |
86: | x x |
87: | x x|
88: p-1 | x|
89: ----------------------------------------
90: */
91: PetscErrorCode DSView_GSVD(DS ds,PetscViewer viewer)
92: {
93: DS_GSVD *ctx = (DS_GSVD*)ds->data;
94: PetscViewerFormat format;
95: PetscInt i,j,r,k=ds->k,n=ds->n,m=ctx->m,p=ctx->p,rowsa,rowsb,colsa,colsb;
96: PetscReal *T,*S,value;
98: PetscFunctionBegin;
99: PetscCall(PetscViewerGetFormat(viewer,&format));
100: if (format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
101: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
102: PetscCall(PetscViewerASCIIPrintf(viewer,"number of columns: %" PetscInt_FMT "\n",m));
103: PetscCall(PetscViewerASCIIPrintf(viewer,"number of rows of B: %" PetscInt_FMT "\n",p));
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
106: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
107: if (ds->compact) {
108: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
109: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
110: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
111: rowsa = n;
112: colsa = ds->extrarow? m+1: m;
113: rowsb = p;
114: colsb = ds->extrarow? m+1: m;
115: if (format == PETSC_VIEWER_ASCII_MATLAB) {
116: PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rowsa,colsa));
117: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",2*ds->n));
118: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
119: for (i=0;i<PetscMin(rowsa,colsa);i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)T[i]));
120: for (i=0;i<k;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,k+1,(double)T[i+ds->ld]));
121: if (n>m) { /* A lower bidiagonal */
122: for (i=k;i<rowsa-1;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+2,i+1,(double)T[i+ds->ld]));
123: } else { /* A (square) upper bidiagonal */
124: for (i=k;i<colsa-1;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+2,(double)T[i+ds->ld]));
125: }
126: PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_T]));
127: PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rowsb,colsb));
128: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",2*ds->n));
129: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
130: for (i=0;i<rowsb;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)S[i]));
131: for (i=0;i<colsb-1;i++) {
132: r = PetscMax(i+2,ds->k+1);
133: PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,r,(double)T[i+2*ds->ld]));
134: }
135: PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_D]));
136: } else {
137: PetscCall(PetscViewerASCIIPrintf(viewer,"Matrix %s =\n",DSMatName[DS_MAT_T]));
138: for (i=0;i<rowsa;i++) {
139: for (j=0;j<colsa;j++) {
140: if (i==j) value = T[i];
141: else if (i<ds->k && j==ds->k) value = T[i+ds->ld];
142: else if (n>m && i==j+1 && i>ds->k) value = T[j+ds->ld];
143: else if (n<=m && i+1==j && i>=ds->k) value = T[i+ds->ld];
144: else value = 0.0;
145: PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
146: }
147: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
148: }
149: PetscCall(PetscViewerASCIIPrintf(viewer,"Matrix %s =\n",DSMatName[DS_MAT_D]));
150: for (i=0;i<rowsb;i++) {
151: for (j=0;j<colsb;j++) {
152: if (i==j) value = S[i];
153: else if (i<ds->k && j==ds->k) value = T[PetscMin(i,j)+2*ds->ld];
154: else if (i+1==j && i>=ds->k) value = T[i+2*ds->ld];
155: else value = 0.0;
156: PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
157: }
158: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
159: }
160: }
161: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
162: PetscCall(PetscViewerFlush(viewer));
163: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
164: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
165: } else {
166: PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
167: PetscCall(DSViewMat(ds,viewer,DS_MAT_B));
168: }
169: if (ds->state>DS_STATE_INTERMEDIATE) {
170: PetscCall(DSViewMat(ds,viewer,DS_MAT_X));
171: PetscCall(DSViewMat(ds,viewer,DS_MAT_U));
172: PetscCall(DSViewMat(ds,viewer,DS_MAT_V));
173: }
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: PetscErrorCode DSVectors_GSVD(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
178: {
179: PetscFunctionBegin;
180: switch (mat) {
181: case DS_MAT_U:
182: case DS_MAT_V:
183: if (rnorm) *rnorm = 0.0;
184: break;
185: case DS_MAT_X:
186: break;
187: default:
188: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
189: }
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: PetscErrorCode DSSort_GSVD(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
194: {
195: DS_GSVD *ctx = (DS_GSVD*)ds->data;
196: PetscInt t,l,ld=ds->ld,i,*perm,*perm2;
197: PetscReal *T=NULL,*D=NULL,*eig;
198: PetscScalar *A=NULL,*B=NULL;
199: PetscBool compact=ds->compact;
201: PetscFunctionBegin;
202: if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
203: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
204: l = ds->l;
205: t = ds->t;
206: perm = ds->perm;
207: PetscCall(PetscMalloc2(t,&eig,t,&perm2));
208: if (compact) {
209: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
210: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
211: for (i=0;i<t;i++) eig[i] = (D[i]==0)?PETSC_INFINITY:T[i]/D[i];
212: } else {
213: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
214: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
215: for (i=0;i<t;i++) eig[i] = (B[i+i*ld]==0)?PETSC_INFINITY:PetscRealPart(A[i+i*ld])/PetscRealPart(B[i*(1+ld)]);
216: }
217: PetscCall(DSSortEigenvaluesReal_Private(ds,eig,perm));
218: PetscCall(PetscArraycpy(perm2,perm,t));
219: for (i=l;i<t;i++) wr[i] = eig[perm[i]];
220: if (compact) {
221: PetscCall(PetscArraycpy(eig,T,t));
222: for (i=l;i<t;i++) T[i] = eig[perm[i]];
223: PetscCall(PetscArraycpy(eig,D,t));
224: for (i=l;i<t;i++) D[i] = eig[perm[i]];
225: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
226: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
227: } else {
228: for (i=l;i<t;i++) eig[i] = PetscRealPart(A[i*(1+ld)]);
229: for (i=l;i<t;i++) A[i*(1+ld)] = eig[perm[i]];
230: for (i=l;i<t;i++) eig[i] = PetscRealPart(B[i*(1+ld)]);
231: for (i=l;i<t;i++) B[i*(1+ld)] = eig[perm[i]];
232: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
233: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
234: }
235: PetscCall(DSPermuteColumns_Private(ds,l,t,ds->n,DS_MAT_U,perm2));
236: PetscCall(PetscArraycpy(perm2,perm,t));
237: PetscCall(DSPermuteColumns_Private(ds,l,t,ctx->m,DS_MAT_X,perm2));
238: PetscCall(DSPermuteColumns_Private(ds,l,t,ctx->p,DS_MAT_V,perm));
239: PetscCall(PetscFree2(eig,perm2));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: PetscErrorCode DSUpdateExtraRow_GSVD(DS ds)
244: {
245: DS_GSVD *ctx = (DS_GSVD*)ds->data;
246: PetscInt i;
247: PetscBLASInt n=0,m=0,ld=0;
248: const PetscScalar *U,*V;
249: PetscReal *T,*e,*f,alpha,beta,betah;
251: PetscFunctionBegin;
252: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
253: PetscCheck(ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-compact storage");
254: PetscCall(PetscBLASIntCast(ds->n,&n));
255: PetscCall(PetscBLASIntCast(ctx->m,&m));
256: PetscCall(PetscBLASIntCast(ds->ld,&ld));
257: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
258: e = T+ld;
259: f = T+2*ld;
260: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_U],&U));
261: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_V],&V));
262: if (n<=m) { /* upper variant, A is square upper bidiagonal */
263: beta = e[m-1]; /* in compact, we assume all entries are zero except the last one */
264: betah = f[m-1];
265: for (i=0;i<m;i++) {
266: e[i] = PetscRealPart(beta*U[m-1+i*ld]);
267: f[i] = PetscRealPart(betah*V[m-1+i*ld]);
268: }
269: } else { /* lower variant, A is (m+1)xm lower bidiagonal */
270: alpha = T[m];
271: betah = f[m-1];
272: for (i=0;i<m;i++) {
273: e[i] = PetscRealPart(alpha*U[m+i*ld]);
274: f[i] = PetscRealPart(betah*V[m-1+i*ld]);
275: }
276: T[m] = PetscRealPart(alpha*U[m+m*ld]);
277: }
278: ds->k = m;
279: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_U],&U));
280: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_V],&V));
281: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: PetscErrorCode DSTruncate_GSVD(DS ds,PetscInt n,PetscBool trim)
286: {
287: DS_GSVD *ctx = (DS_GSVD*)ds->data;
288: PetscScalar *U;
289: PetscReal *T;
290: PetscInt i,m=ctx->m,ld=ds->ld;
291: PetscBool lower=(ds->n>ctx->m)?PETSC_TRUE:PETSC_FALSE;
293: PetscFunctionBegin;
294: PetscCheck(ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-compact storage");
295: if (trim) {
296: ds->l = 0;
297: ds->k = 0;
298: ds->n = lower? n+1: n;
299: ctx->m = n;
300: ctx->p = n;
301: ds->t = ds->n; /* truncated length equal to the new dimension */
302: ctx->tm = ctx->m; /* must also keep the previous dimension of X */
303: ctx->tp = ctx->p; /* must also keep the previous dimension of V */
304: } else {
305: if (lower) {
306: /* move value of diagonal element of arrow (alpha) */
307: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
308: T[n] = T[m];
309: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
310: /* copy last column of U so that it updates the next initial vector of U1 */
311: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
312: for (i=0;i<=m;i++) U[i+n*ld] = U[i+m*ld];
313: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
314: }
315: ds->k = (ds->extrarow)? n: 0;
316: ds->t = ds->n; /* truncated length equal to previous dimension */
317: ctx->tm = ctx->m; /* must also keep the previous dimension of X */
318: ctx->tp = ctx->p; /* must also keep the previous dimension of V */
319: ds->n = lower? n+1: n;
320: ctx->m = n;
321: ctx->p = n;
322: }
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: static PetscErrorCode DSSwitchFormat_GSVD(DS ds)
327: {
328: DS_GSVD *ctx = (DS_GSVD*)ds->data;
329: PetscReal *T,*D;
330: PetscScalar *A,*B;
331: PetscInt i,n=ds->n,k=ds->k,ld=ds->ld,m=ctx->m;
333: PetscFunctionBegin;
334: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
335: /* switch from compact (arrow) to dense storage */
336: /* bidiagonal associated to B is stored in D and T+2*ld */
337: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_A],&A));
338: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_B],&B));
339: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
340: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
341: PetscCall(PetscArrayzero(A,ld*ld));
342: PetscCall(PetscArrayzero(B,ld*ld));
343: for (i=0;i<k;i++) {
344: A[i+i*ld] = T[i];
345: A[i+k*ld] = T[i+ld];
346: B[i+i*ld] = D[i];
347: B[i+k*ld] = T[i+2*ld];
348: }
349: /* B is upper bidiagonal */
350: B[k+k*ld] = D[k];
351: for (i=k+1;i<m;i++) {
352: B[i+i*ld] = D[i];
353: B[i-1+i*ld] = T[i-1+2*ld];
354: }
355: /* A can be upper (square) or lower bidiagonal */
356: for (i=k;i<m;i++) A[i+i*ld] = T[i];
357: if (n>m) for (i=k;i<m;i++) A[i+1+i*ld] = T[i+ld];
358: else for (i=k+1;i<m;i++) A[i-1+i*ld] = T[i-1+ld];
359: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_A],&A));
360: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_B],&B));
361: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
362: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: /*
367: Compact format is used when [A;B] has orthonormal columns.
368: In this case R=I and the GSVD of (A,B) is the CS decomposition
369: */
370: PetscErrorCode DSSolve_GSVD(DS ds,PetscScalar *wr,PetscScalar *wi)
371: {
372: DS_GSVD *ctx = (DS_GSVD*)ds->data;
373: PetscInt i,j;
374: PetscBLASInt n1,m1,info,lc = 0,n = 0,m = 0,p = 0,p1,l,k,q,ld,off,lwork,r;
375: PetscScalar *A,*B,*X,*U,*V,sone=1.0,smone=-1.0;
376: PetscReal *alpha,*beta,*T,*D;
377: #if !defined(SLEPC_MISSING_LAPACK_GGSVD3)
378: PetscScalar a,dummy;
379: PetscReal rdummy;
380: PetscBLASInt idummy;
381: #endif
383: PetscFunctionBegin;
384: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
385: PetscCall(PetscBLASIntCast(ds->n,&m));
386: PetscCall(PetscBLASIntCast(ctx->m,&n));
387: PetscCall(PetscBLASIntCast(ctx->p,&p));
388: PetscCall(PetscBLASIntCast(ds->l,&lc));
389: PetscCheck(ds->compact || lc==0,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DSGSVD with non-compact format does not support locking");
390: /* In compact storage B is always nxn and A can be either nxn or (n+1)xn */
391: PetscCheck(!ds->compact || (p==n && (m==p || m==p+1)),PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Dimensions not supported in compact format");
392: PetscCall(PetscBLASIntCast(ds->ld,&ld));
393: n1 = n-lc; /* n1 = size of leading block, excl. locked + size of trailing block */
394: m1 = m-lc;
395: p1 = p-lc;
396: off = lc+lc*ld;
397: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
398: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
399: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
400: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
401: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
402: PetscCall(PetscArrayzero(X,ld*ld));
403: for (i=0;i<lc;i++) X[i+i*ld] = 1.0;
404: PetscCall(PetscArrayzero(U,ld*ld));
405: for (i=0;i<lc;i++) U[i+i*ld] = 1.0;
406: PetscCall(PetscArrayzero(V,ld*ld));
407: for (i=0;i<lc;i++) V[i+i*ld] = 1.0;
408: if (ds->compact) PetscCall(DSSwitchFormat_GSVD(ds));
410: #if !defined(SLEPC_MISSING_LAPACK_GGSVD3)
411: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
412: /* workspace query and memory allocation */
413: lwork = -1;
414: #if !defined (PETSC_USE_COMPLEX)
415: PetscCallBLAS("LAPACKggsvd3",LAPACKggsvd3_("U","V","Q",&m1,&n1,&p1,&k,&l,&dummy,&ld,&dummy,&ld,&rdummy,&rdummy,&dummy,&ld,&dummy,&ld,&dummy,&ld,&a,&lwork,&idummy,&info));
416: PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
417: #else
418: PetscCallBLAS("LAPACKggsvd3",LAPACKggsvd3_("U","V","Q",&m1,&n1,&p1,&k,&l,&dummy,&ld,&dummy,&ld,&rdummy,&rdummy,&dummy,&ld,&dummy,&ld,&dummy,&ld,&a,&lwork,&rdummy,&idummy,&info));
419: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
420: #endif
422: #if !defined (PETSC_USE_COMPLEX)
423: PetscCall(DSAllocateWork_Private(ds,lwork,2*ds->ld,ds->ld));
424: alpha = ds->rwork;
425: beta = ds->rwork+ds->ld;
426: PetscCallBLAS("LAPACKggsvd3",LAPACKggsvd3_("U","V","Q",&m1,&n1,&p1,&k,&l,A+off,&ld,B+off,&ld,alpha,beta,U+off,&ld,V+off,&ld,X+off,&ld,ds->work,&lwork,ds->iwork,&info));
427: #else
428: PetscCall(DSAllocateWork_Private(ds,lwork,4*ds->ld,ds->ld));
429: alpha = ds->rwork+2*ds->ld;
430: beta = ds->rwork+3*ds->ld;
431: PetscCallBLAS("LAPACKggsvd3",LAPACKggsvd3_("U","V","Q",&m1,&n1,&p1,&k,&l,A+off,&ld,B+off,&ld,alpha,beta,U+off,&ld,V+off,&ld,X+off,&ld,ds->work,&lwork,ds->rwork,ds->iwork,&info));
432: #endif
433: PetscCall(PetscFPTrapPop());
434: SlepcCheckLapackInfo("ggsvd3",info);
436: #else /* defined(SLEPC_MISSING_LAPACK_GGSVD3) */
438: lwork = PetscMax(PetscMax(3*n,m),p)+n;
439: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
440: #if !defined (PETSC_USE_COMPLEX)
441: PetscCall(DSAllocateWork_Private(ds,lwork,2*ds->ld,ds->ld));
442: alpha = ds->rwork;
443: beta = ds->rwork+ds->ld;
444: PetscCallBLAS("LAPACKggsvd",LAPACKggsvd_("U","V","Q",&m1,&n1,&p1,&k,&l,A+off,&ld,B+off,&ld,alpha,beta,U+off,&ld,V+off,&ld,X+off,&ld,ds->work,ds->iwork,&info));
445: #else
446: PetscCall(DSAllocateWork_Private(ds,lwork,4*ds->ld,ds->ld));
447: alpha = ds->rwork+2*ds->ld;
448: beta = ds->rwork+3*ds->ld;
449: PetscCallBLAS("LAPACKggsvd",LAPACKggsvd_("U","V","Q",&m1,&n1,&p1,&k,&l,A+off,&ld,B+off,&ld,alpha,beta,U+off,&ld,V+off,&ld,X+off,&ld,ds->work,ds->rwork,ds->iwork,&info));
450: #endif
451: PetscCall(PetscFPTrapPop());
452: SlepcCheckLapackInfo("ggsvd",info);
454: #endif
456: PetscCheck(k+l>=n1,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"The rank deficient case not supported yet");
457: if (ds->compact) {
458: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
459: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
460: /* R is the identity matrix (except the sign) */
461: for (i=lc;i<n;i++) {
462: if (PetscRealPart(A[i+i*ld])<0.0) { /* scale column i */
463: for (j=lc;j<n;j++) X[j+i*ld] = -X[j+i*ld];
464: }
465: }
466: PetscCall(PetscArrayzero(T+ld,m-1));
467: PetscCall(PetscArrayzero(T+2*ld,n-1));
468: for (i=lc;i<n;i++) {
469: T[i] = alpha[i-lc];
470: D[i] = beta[i-lc];
471: if (D[i]==0.0) wr[i] = PETSC_INFINITY;
472: else wr[i] = T[i]/D[i];
473: }
474: ds->t = n;
475: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
476: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
477: } else {
478: /* X = X*inv(R) */
479: q = PetscMin(m,n);
480: PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n,&q,&sone,A,&ld,X,&ld));
481: if (m<n) {
482: r = n-m;
483: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&r,&m,&sone,X,&ld,A,&ld,&smone,X+m*ld,&ld));
484: PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n,&r,&sone,B+m*ld,&ld,X+m*ld,&ld));
485: }
486: if (k>0) {
487: for (i=k;i<PetscMin(m,k+l);i++) {
488: PetscCall(PetscArraycpy(X+(i-k)*ld,X+i*ld,ld));
489: PetscCall(PetscArraycpy(U+(i-k)*ld,U+i*ld,ld));
490: }
491: }
492: /* singular values */
493: PetscCall(PetscArrayzero(A,ld*ld));
494: PetscCall(PetscArrayzero(B,ld*ld));
495: for (j=k;j<PetscMin(m,k+l);j++) {
496: A[(j-k)*(1+ld)] = alpha[j];
497: B[(j-k)*(1+ld)] = beta[j];
498: wr[j-k] = alpha[j]/beta[j];
499: }
500: ds->t = PetscMin(m,k+l)-k; /* set number of computed values */
501: }
502: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
503: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
504: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
505: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
506: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }
510: PetscErrorCode DSCond_GSVD(DS ds,PetscReal *cond)
511: {
512: DS_GSVD *ctx = (DS_GSVD*)ds->data;
513: PetscBLASInt lwork,lrwork=0,info,m,n,p,ld;
514: PetscScalar *A,*work;
515: const PetscScalar *M;
516: PetscReal *sigma,conda,condb;
517: #if defined(PETSC_USE_COMPLEX)
518: PetscReal *rwork;
519: #endif
521: PetscFunctionBegin;
522: PetscCall(PetscBLASIntCast(ds->n,&m));
523: PetscCall(PetscBLASIntCast(ctx->m,&n));
524: PetscCall(PetscBLASIntCast(ctx->p,&p));
525: PetscCall(PetscBLASIntCast(ds->ld,&ld));
526: lwork = 5*n;
527: #if defined(PETSC_USE_COMPLEX)
528: lrwork = 5*n;
529: #endif
530: PetscCall(DSAllocateWork_Private(ds,ld*n+lwork,n+lrwork,0));
531: A = ds->work;
532: work = ds->work+ld*n;
533: sigma = ds->rwork;
534: #if defined(PETSC_USE_COMPLEX)
535: rwork = ds->rwork+n;
536: #endif
537: if (ds->compact) PetscCall(DSSwitchFormat_GSVD(ds));
539: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
540: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&M));
541: PetscCall(PetscArraycpy(A,M,ld*n));
542: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&M));
543: #if defined(PETSC_USE_COMPLEX)
544: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,rwork,&info));
545: #else
546: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,&info));
547: #endif
548: SlepcCheckLapackInfo("gesvd",info);
549: conda = sigma[0]/sigma[PetscMin(m,n)-1];
551: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&M));
552: PetscCall(PetscArraycpy(A,M,ld*n));
553: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&M));
554: #if defined(PETSC_USE_COMPLEX)
555: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&p,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,rwork,&info));
556: #else
557: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&p,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,&info));
558: #endif
559: SlepcCheckLapackInfo("gesvd",info);
560: condb = sigma[0]/sigma[PetscMin(p,n)-1];
561: PetscCall(PetscFPTrapPop());
563: *cond = PetscMax(conda,condb);
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: #if !defined(PETSC_HAVE_MPIUNI)
568: PetscErrorCode DSSynchronize_GSVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
569: {
570: DS_GSVD *ctx = (DS_GSVD*)ds->data;
571: PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
572: PetscMPIInt m=ctx->m,rank,off=0,size,n,ldn,ld3;
573: PetscScalar *A,*U,*V,*X;
574: PetscReal *T;
576: PetscFunctionBegin;
577: if (ds->compact) kr = 3*ld;
578: else k = 2*(m-l)*ld;
579: if (ds->state>DS_STATE_RAW) k += 3*(m-l)*ld;
580: if (eigr) k += m-l;
581: PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
582: PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
583: PetscCall(PetscMPIIntCast(m-l,&n));
584: PetscCall(PetscMPIIntCast(ld*(m-l),&ldn));
585: PetscCall(PetscMPIIntCast(3*ld,&ld3));
586: if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
587: else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
588: if (ds->state>DS_STATE_RAW) {
589: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
590: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
591: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
592: }
593: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
594: if (!rank) {
595: if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
596: else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
597: if (ds->state>DS_STATE_RAW) {
598: PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
599: PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
600: PetscCallMPI(MPI_Pack(X+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
601: }
602: if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
603: }
604: PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
605: if (rank) {
606: if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
607: else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
608: if (ds->state>DS_STATE_RAW) {
609: PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
610: PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
611: PetscCallMPI(MPI_Unpack(ds->work,size,&off,X+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
612: }
613: if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
614: }
615: if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
616: else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
617: if (ds->state>DS_STATE_RAW) {
618: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
619: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
620: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
621: }
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
624: #endif
626: PetscErrorCode DSMatGetSize_GSVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
627: {
628: DS_GSVD *ctx = (DS_GSVD*)ds->data;
630: PetscFunctionBegin;
631: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
632: switch (t) {
633: case DS_MAT_A:
634: *rows = ds->n;
635: *cols = ds->extrarow? ctx->m+1: ctx->m;
636: break;
637: case DS_MAT_B:
638: *rows = ctx->p;
639: *cols = ds->extrarow? ctx->m+1: ctx->m;
640: break;
641: case DS_MAT_T:
642: *rows = ds->n;
643: *cols = PetscDefined(USE_COMPLEX)? 2: 3;
644: break;
645: case DS_MAT_D:
646: *rows = ctx->p;
647: *cols = 1;
648: break;
649: case DS_MAT_U:
650: *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
651: *cols = ds->n;
652: break;
653: case DS_MAT_V:
654: *rows = ds->state==DS_STATE_TRUNCATED? ctx->tp: ctx->p;
655: *cols = ctx->p;
656: break;
657: case DS_MAT_X:
658: *rows = ds->state==DS_STATE_TRUNCATED? ctx->tm: ctx->m;
659: *cols = ctx->m;
660: break;
661: default:
662: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
663: }
664: PetscFunctionReturn(PETSC_SUCCESS);
665: }
667: static PetscErrorCode DSGSVDSetDimensions_GSVD(DS ds,PetscInt m,PetscInt p)
668: {
669: DS_GSVD *ctx = (DS_GSVD*)ds->data;
671: PetscFunctionBegin;
672: DSCheckAlloc(ds,1);
673: if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
674: ctx->m = ds->ld;
675: } else {
676: PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
677: ctx->m = m;
678: }
679: if (p==PETSC_DECIDE || p==PETSC_DEFAULT) {
680: ctx->p = ds->n;
681: } else {
682: PetscCheck(p>0 && p<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of p. Must be between 1 and ld");
683: ctx->p = p;
684: }
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
688: /*@
689: DSGSVDSetDimensions - Sets the number of columns and rows for a DSGSVD.
691: Logically Collective
693: Input Parameters:
694: + ds - the direct solver context
695: . m - the number of columns
696: - p - the number of rows for the second matrix (B)
698: Notes:
699: This call is complementary to DSSetDimensions(), to provide two dimensions
700: that are specific to this DS type. The number of rows for the first matrix (A)
701: is set by DSSetDimensions().
703: Level: intermediate
705: .seealso: DSGSVDGetDimensions(), DSSetDimensions()
706: @*/
707: PetscErrorCode DSGSVDSetDimensions(DS ds,PetscInt m,PetscInt p)
708: {
709: PetscFunctionBegin;
713: PetscTryMethod(ds,"DSGSVDSetDimensions_C",(DS,PetscInt,PetscInt),(ds,m,p));
714: PetscFunctionReturn(PETSC_SUCCESS);
715: }
717: static PetscErrorCode DSGSVDGetDimensions_GSVD(DS ds,PetscInt *m,PetscInt *p)
718: {
719: DS_GSVD *ctx = (DS_GSVD*)ds->data;
721: PetscFunctionBegin;
722: if (m) *m = ctx->m;
723: if (p) *p = ctx->p;
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: /*@
728: DSGSVDGetDimensions - Returns the number of columns and rows for a DSGSVD.
730: Not Collective
732: Input Parameter:
733: . ds - the direct solver context
735: Output Parameters:
736: + m - the number of columns
737: - p - the number of rows for the second matrix (B)
739: Level: intermediate
741: .seealso: DSGSVDSetDimensions()
742: @*/
743: PetscErrorCode DSGSVDGetDimensions(DS ds,PetscInt *m,PetscInt *p)
744: {
745: PetscFunctionBegin;
747: PetscUseMethod(ds,"DSGSVDGetDimensions_C",(DS,PetscInt*,PetscInt*),(ds,m,p));
748: PetscFunctionReturn(PETSC_SUCCESS);
749: }
751: PetscErrorCode DSDestroy_GSVD(DS ds)
752: {
753: PetscFunctionBegin;
754: PetscCall(PetscFree(ds->data));
755: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDSetDimensions_C",NULL));
756: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDGetDimensions_C",NULL));
757: PetscFunctionReturn(PETSC_SUCCESS);
758: }
760: /*MC
761: DSGSVD - Dense Generalized Singular Value Decomposition.
763: Level: beginner
765: Notes:
766: The problem is expressed as A*X = U*C, B*X = V*S, where A and B are
767: matrices with the same number of columns, m, U and V are orthogonal
768: (unitary), and X is an mxm invertible matrix. The DS object does not
769: expose matrices C and S, instead the singular values sigma, which are
770: the ratios c_i/s_i, are returned in the arguments of DSSolve().
771: Note that the number of columns of the returned X, U, V may be smaller
772: in the case that some c_i or s_i are zero.
774: The number of rows of A (and U) is the value n passed with DSSetDimensions().
775: The number of columns m and the number of rows of B (and V) must be
776: set via DSGSVDSetDimensions().
778: Internally, LAPACK's representation is used, U'*A*Q = C*[0 R], V'*B*Q = S*[0 R],
779: where X = Q*inv(R) is computed at the end of DSSolve().
781: If the compact storage format is selected, then a simplified problem is
782: solved, where A and B are bidiagonal (possibly with an arrow), and [A;B]
783: is assumed to have orthonormal columns. We consider two cases: (1) A and B
784: are square mxm upper bidiagonal, and (2) A is lower bidiagonal with m+1
785: rows and B is square upper bidiagonal. In these cases, R=I so it
786: corresponds to the CS decomposition. The first matrix is stored in two
787: diagonals of DS_MAT_T, while the second matrix is stored in DS_MAT_D
788: and the remaining diagonal of DS_MAT_T.
790: Allowed arguments of DSVectors() are DS_MAT_U, DS_MAT_V and DS_MAT_X.
792: Used DS matrices:
793: + DS_MAT_A - first problem matrix
794: . DS_MAT_B - second problem matrix
795: . DS_MAT_T - first upper bidiagonal matrix (if compact storage is selected)
796: - DS_MAT_D - second upper bidiagonal matrix (if compact storage is selected)
798: Implemented methods:
799: . 0 - Lapack (_ggsvd3 if available, or _ggsvd)
801: .seealso: DSCreate(), DSSetType(), DSType, DSGSVDSetDimensions()
802: M*/
803: SLEPC_EXTERN PetscErrorCode DSCreate_GSVD(DS ds)
804: {
805: DS_GSVD *ctx;
807: PetscFunctionBegin;
808: PetscCall(PetscNew(&ctx));
809: ds->data = (void*)ctx;
811: ds->ops->allocate = DSAllocate_GSVD;
812: ds->ops->view = DSView_GSVD;
813: ds->ops->vectors = DSVectors_GSVD;
814: ds->ops->sort = DSSort_GSVD;
815: ds->ops->solve[0] = DSSolve_GSVD;
816: #if !defined(PETSC_HAVE_MPIUNI)
817: ds->ops->synchronize = DSSynchronize_GSVD;
818: #endif
819: ds->ops->truncate = DSTruncate_GSVD;
820: ds->ops->update = DSUpdateExtraRow_GSVD;
821: ds->ops->cond = DSCond_GSVD;
822: ds->ops->matgetsize = DSMatGetSize_GSVD;
823: ds->ops->destroy = DSDestroy_GSVD;
824: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDSetDimensions_C",DSGSVDSetDimensions_GSVD));
825: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDGetDimensions_C",DSGSVDGetDimensions_GSVD));
826: PetscFunctionReturn(PETSC_SUCCESS);
827: }