Actual source code: dsgsvd.c
slepc-3.21.1 2024-04-26
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: static 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: static 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: static 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: static 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: static 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: static 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: static 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: /* workspace query and memory allocation */
412: lwork = -1;
413: #if !defined (PETSC_USE_COMPLEX)
414: 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));
415: PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
416: #else
417: 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));
418: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
419: #endif
421: #if !defined (PETSC_USE_COMPLEX)
422: PetscCall(DSAllocateWork_Private(ds,lwork,2*ds->ld,ds->ld));
423: alpha = ds->rwork;
424: beta = ds->rwork+ds->ld;
425: 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));
426: #else
427: PetscCall(DSAllocateWork_Private(ds,lwork,4*ds->ld,ds->ld));
428: alpha = ds->rwork+2*ds->ld;
429: beta = ds->rwork+3*ds->ld;
430: 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));
431: #endif
432: SlepcCheckLapackInfo("ggsvd3",info);
434: #else /* defined(SLEPC_MISSING_LAPACK_GGSVD3) */
436: lwork = PetscMax(PetscMax(3*n,m),p)+n;
437: #if !defined (PETSC_USE_COMPLEX)
438: PetscCall(DSAllocateWork_Private(ds,lwork,2*ds->ld,ds->ld));
439: alpha = ds->rwork;
440: beta = ds->rwork+ds->ld;
441: 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));
442: #else
443: PetscCall(DSAllocateWork_Private(ds,lwork,4*ds->ld,ds->ld));
444: alpha = ds->rwork+2*ds->ld;
445: beta = ds->rwork+3*ds->ld;
446: 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));
447: #endif
448: SlepcCheckLapackInfo("ggsvd",info);
450: #endif
452: PetscCheck(k+l>=n1,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"The rank deficient case not supported yet");
453: if (ds->compact) {
454: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
455: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
456: /* R is the identity matrix (except the sign) */
457: for (i=lc;i<n;i++) {
458: if (PetscRealPart(A[i+i*ld])<0.0) { /* scale column i */
459: for (j=lc;j<n;j++) X[j+i*ld] = -X[j+i*ld];
460: }
461: }
462: PetscCall(PetscArrayzero(T+ld,m-1));
463: PetscCall(PetscArrayzero(T+2*ld,n-1));
464: for (i=lc;i<n;i++) {
465: T[i] = alpha[i-lc];
466: D[i] = beta[i-lc];
467: if (D[i]==0.0) wr[i] = PETSC_INFINITY;
468: else wr[i] = T[i]/D[i];
469: }
470: ds->t = n;
471: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
472: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
473: } else {
474: /* X = X*inv(R) */
475: q = PetscMin(m,n);
476: PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n,&q,&sone,A,&ld,X,&ld));
477: if (m<n) {
478: r = n-m;
479: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&r,&m,&sone,X,&ld,A,&ld,&smone,X+m*ld,&ld));
480: PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n,&r,&sone,B+m*ld,&ld,X+m*ld,&ld));
481: }
482: if (k>0) {
483: for (i=k;i<PetscMin(m,k+l);i++) {
484: PetscCall(PetscArraycpy(X+(i-k)*ld,X+i*ld,ld));
485: PetscCall(PetscArraycpy(U+(i-k)*ld,U+i*ld,ld));
486: }
487: }
488: /* singular values */
489: PetscCall(PetscArrayzero(A,ld*ld));
490: PetscCall(PetscArrayzero(B,ld*ld));
491: for (j=k;j<PetscMin(m,k+l);j++) {
492: A[(j-k)*(1+ld)] = alpha[j];
493: B[(j-k)*(1+ld)] = beta[j];
494: wr[j-k] = alpha[j]/beta[j];
495: }
496: ds->t = PetscMin(m,k+l)-k; /* set number of computed values */
497: }
498: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
499: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
500: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
501: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
502: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: static PetscErrorCode DSCond_GSVD(DS ds,PetscReal *cond)
507: {
508: DS_GSVD *ctx = (DS_GSVD*)ds->data;
509: PetscBLASInt lwork,lrwork=0,info,m,n,p,ld;
510: PetscScalar *A,*work;
511: const PetscScalar *M;
512: PetscReal *sigma,conda,condb;
513: #if defined(PETSC_USE_COMPLEX)
514: PetscReal *rwork;
515: #endif
517: PetscFunctionBegin;
518: PetscCall(PetscBLASIntCast(ds->n,&m));
519: PetscCall(PetscBLASIntCast(ctx->m,&n));
520: PetscCall(PetscBLASIntCast(ctx->p,&p));
521: PetscCall(PetscBLASIntCast(ds->ld,&ld));
522: lwork = 5*n;
523: #if defined(PETSC_USE_COMPLEX)
524: lrwork = 5*n;
525: #endif
526: PetscCall(DSAllocateWork_Private(ds,ld*n+lwork,n+lrwork,0));
527: A = ds->work;
528: work = ds->work+ld*n;
529: sigma = ds->rwork;
530: #if defined(PETSC_USE_COMPLEX)
531: rwork = ds->rwork+n;
532: #endif
533: if (ds->compact) PetscCall(DSSwitchFormat_GSVD(ds));
535: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&M));
536: PetscCall(PetscArraycpy(A,M,ld*n));
537: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&M));
538: #if defined(PETSC_USE_COMPLEX)
539: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,rwork,&info));
540: #else
541: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,&info));
542: #endif
543: SlepcCheckLapackInfo("gesvd",info);
544: conda = sigma[0]/sigma[PetscMin(m,n)-1];
546: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&M));
547: PetscCall(PetscArraycpy(A,M,ld*n));
548: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&M));
549: #if defined(PETSC_USE_COMPLEX)
550: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&p,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,rwork,&info));
551: #else
552: PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&p,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,&info));
553: #endif
554: SlepcCheckLapackInfo("gesvd",info);
555: condb = sigma[0]/sigma[PetscMin(p,n)-1];
557: *cond = PetscMax(conda,condb);
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: #if !defined(PETSC_HAVE_MPIUNI)
562: static PetscErrorCode DSSynchronize_GSVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
563: {
564: DS_GSVD *ctx = (DS_GSVD*)ds->data;
565: PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
566: PetscMPIInt m=ctx->m,rank,off=0,size,n,ldn,ld3;
567: PetscScalar *A,*U,*V,*X;
568: PetscReal *T;
570: PetscFunctionBegin;
571: if (ds->compact) kr = 3*ld;
572: else k = 2*(m-l)*ld;
573: if (ds->state>DS_STATE_RAW) k += 3*(m-l)*ld;
574: if (eigr) k += m-l;
575: PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
576: PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
577: PetscCall(PetscMPIIntCast(m-l,&n));
578: PetscCall(PetscMPIIntCast(ld*(m-l),&ldn));
579: PetscCall(PetscMPIIntCast(3*ld,&ld3));
580: if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
581: else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
582: if (ds->state>DS_STATE_RAW) {
583: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
584: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
585: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
586: }
587: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
588: if (!rank) {
589: if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
590: else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
591: if (ds->state>DS_STATE_RAW) {
592: PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
593: PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
594: PetscCallMPI(MPI_Pack(X+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
595: }
596: if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
597: }
598: PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
599: if (rank) {
600: if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
601: else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
602: if (ds->state>DS_STATE_RAW) {
603: PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
604: PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
605: PetscCallMPI(MPI_Unpack(ds->work,size,&off,X+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
606: }
607: if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
608: }
609: if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
610: else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
611: if (ds->state>DS_STATE_RAW) {
612: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
613: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
614: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
615: }
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
618: #endif
620: static PetscErrorCode DSMatGetSize_GSVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
621: {
622: DS_GSVD *ctx = (DS_GSVD*)ds->data;
624: PetscFunctionBegin;
625: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
626: switch (t) {
627: case DS_MAT_A:
628: *rows = ds->n;
629: *cols = ds->extrarow? ctx->m+1: ctx->m;
630: break;
631: case DS_MAT_B:
632: *rows = ctx->p;
633: *cols = ds->extrarow? ctx->m+1: ctx->m;
634: break;
635: case DS_MAT_T:
636: *rows = ds->n;
637: *cols = PetscDefined(USE_COMPLEX)? 2: 3;
638: break;
639: case DS_MAT_D:
640: *rows = ctx->p;
641: *cols = 1;
642: break;
643: case DS_MAT_U:
644: *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
645: *cols = ds->n;
646: break;
647: case DS_MAT_V:
648: *rows = ds->state==DS_STATE_TRUNCATED? ctx->tp: ctx->p;
649: *cols = ctx->p;
650: break;
651: case DS_MAT_X:
652: *rows = ds->state==DS_STATE_TRUNCATED? ctx->tm: ctx->m;
653: *cols = ctx->m;
654: break;
655: default:
656: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
657: }
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }
661: static PetscErrorCode DSGSVDSetDimensions_GSVD(DS ds,PetscInt m,PetscInt p)
662: {
663: DS_GSVD *ctx = (DS_GSVD*)ds->data;
665: PetscFunctionBegin;
666: DSCheckAlloc(ds,1);
667: if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
668: ctx->m = ds->ld;
669: } else {
670: PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
671: ctx->m = m;
672: }
673: if (p==PETSC_DECIDE || p==PETSC_DEFAULT) {
674: ctx->p = ds->n;
675: } else {
676: PetscCheck(p>0 && p<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of p. Must be between 1 and ld");
677: ctx->p = p;
678: }
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: /*@
683: DSGSVDSetDimensions - Sets the number of columns and rows for a DSGSVD.
685: Logically Collective
687: Input Parameters:
688: + ds - the direct solver context
689: . m - the number of columns
690: - p - the number of rows for the second matrix (B)
692: Notes:
693: This call is complementary to DSSetDimensions(), to provide two dimensions
694: that are specific to this DS type. The number of rows for the first matrix (A)
695: is set by DSSetDimensions().
697: Level: intermediate
699: .seealso: DSGSVDGetDimensions(), DSSetDimensions()
700: @*/
701: PetscErrorCode DSGSVDSetDimensions(DS ds,PetscInt m,PetscInt p)
702: {
703: PetscFunctionBegin;
707: PetscTryMethod(ds,"DSGSVDSetDimensions_C",(DS,PetscInt,PetscInt),(ds,m,p));
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: static PetscErrorCode DSGSVDGetDimensions_GSVD(DS ds,PetscInt *m,PetscInt *p)
712: {
713: DS_GSVD *ctx = (DS_GSVD*)ds->data;
715: PetscFunctionBegin;
716: if (m) *m = ctx->m;
717: if (p) *p = ctx->p;
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: /*@
722: DSGSVDGetDimensions - Returns the number of columns and rows for a DSGSVD.
724: Not Collective
726: Input Parameter:
727: . ds - the direct solver context
729: Output Parameters:
730: + m - the number of columns
731: - p - the number of rows for the second matrix (B)
733: Level: intermediate
735: .seealso: DSGSVDSetDimensions()
736: @*/
737: PetscErrorCode DSGSVDGetDimensions(DS ds,PetscInt *m,PetscInt *p)
738: {
739: PetscFunctionBegin;
741: PetscUseMethod(ds,"DSGSVDGetDimensions_C",(DS,PetscInt*,PetscInt*),(ds,m,p));
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }
745: static PetscErrorCode DSDestroy_GSVD(DS ds)
746: {
747: PetscFunctionBegin;
748: PetscCall(PetscFree(ds->data));
749: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDSetDimensions_C",NULL));
750: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDGetDimensions_C",NULL));
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: /*MC
755: DSGSVD - Dense Generalized Singular Value Decomposition.
757: Level: beginner
759: Notes:
760: The problem is expressed as A*X = U*C, B*X = V*S, where A and B are
761: matrices with the same number of columns, m, U and V are orthogonal
762: (unitary), and X is an mxm invertible matrix. The DS object does not
763: expose matrices C and S, instead the singular values sigma, which are
764: the ratios c_i/s_i, are returned in the arguments of DSSolve().
765: Note that the number of columns of the returned X, U, V may be smaller
766: in the case that some c_i or s_i are zero.
768: The number of rows of A (and U) is the value n passed with DSSetDimensions().
769: The number of columns m and the number of rows of B (and V) must be
770: set via DSGSVDSetDimensions().
772: Internally, LAPACK's representation is used, U'*A*Q = C*[0 R], V'*B*Q = S*[0 R],
773: where X = Q*inv(R) is computed at the end of DSSolve().
775: If the compact storage format is selected, then a simplified problem is
776: solved, where A and B are bidiagonal (possibly with an arrow), and [A;B]
777: is assumed to have orthonormal columns. We consider two cases: (1) A and B
778: are square mxm upper bidiagonal, and (2) A is lower bidiagonal with m+1
779: rows and B is square upper bidiagonal. In these cases, R=I so it
780: corresponds to the CS decomposition. The first matrix is stored in two
781: diagonals of DS_MAT_T, while the second matrix is stored in DS_MAT_D
782: and the remaining diagonal of DS_MAT_T.
784: Allowed arguments of DSVectors() are DS_MAT_U, DS_MAT_V and DS_MAT_X.
786: Used DS matrices:
787: + DS_MAT_A - first problem matrix
788: . DS_MAT_B - second problem matrix
789: . DS_MAT_T - first upper bidiagonal matrix (if compact storage is selected)
790: - DS_MAT_D - second upper bidiagonal matrix (if compact storage is selected)
792: Implemented methods:
793: . 0 - Lapack (_ggsvd3 if available, or _ggsvd)
795: .seealso: DSCreate(), DSSetType(), DSType, DSGSVDSetDimensions()
796: M*/
797: SLEPC_EXTERN PetscErrorCode DSCreate_GSVD(DS ds)
798: {
799: DS_GSVD *ctx;
801: PetscFunctionBegin;
802: PetscCall(PetscNew(&ctx));
803: ds->data = (void*)ctx;
805: ds->ops->allocate = DSAllocate_GSVD;
806: ds->ops->view = DSView_GSVD;
807: ds->ops->vectors = DSVectors_GSVD;
808: ds->ops->sort = DSSort_GSVD;
809: ds->ops->solve[0] = DSSolve_GSVD;
810: #if !defined(PETSC_HAVE_MPIUNI)
811: ds->ops->synchronize = DSSynchronize_GSVD;
812: #endif
813: ds->ops->truncate = DSTruncate_GSVD;
814: ds->ops->update = DSUpdateExtraRow_GSVD;
815: ds->ops->cond = DSCond_GSVD;
816: ds->ops->matgetsize = DSMatGetSize_GSVD;
817: ds->ops->destroy = DSDestroy_GSVD;
818: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDSetDimensions_C",DSGSVDSetDimensions_GSVD));
819: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDGetDimensions_C",DSGSVDGetDimensions_GSVD));
820: PetscFunctionReturn(PETSC_SUCCESS);
821: }