Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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> /*I "slepcds.h" I*/
11 : #include <slepcblaslapack.h>
12 :
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;
19 :
20 72 : static PetscErrorCode DSAllocate_GSVD(DS ds,PetscInt ld)
21 : {
22 72 : PetscFunctionBegin;
23 72 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
24 72 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
25 72 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
26 72 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
27 72 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_V));
28 72 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
29 72 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_D));
30 72 : PetscCall(PetscFree(ds->perm));
31 72 : PetscCall(PetscMalloc1(ld,&ds->perm));
32 72 : PetscFunctionReturn(PETSC_SUCCESS);
33 : }
34 :
35 : /*
36 : In compact form, A is either in form (a) or (b):
37 :
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 : ----------------------------------------- -----------------------------------------
64 :
65 : and B is square bidiagonal with upper arrow (p=m)
66 :
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 30 : static PetscErrorCode DSView_GSVD(DS ds,PetscViewer viewer)
92 : {
93 30 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
94 30 : PetscViewerFormat format;
95 30 : PetscInt i,j,r,k=ds->k,n=ds->n,m=ctx->m,p=ctx->p,rowsa,rowsb,colsa,colsb;
96 30 : PetscReal *T,*S,value;
97 :
98 30 : PetscFunctionBegin;
99 30 : PetscCall(PetscViewerGetFormat(viewer,&format));
100 30 : if (format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
101 30 : if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
102 17 : PetscCall(PetscViewerASCIIPrintf(viewer,"number of columns: %" PetscInt_FMT "\n",m));
103 17 : PetscCall(PetscViewerASCIIPrintf(viewer,"number of rows of B: %" PetscInt_FMT "\n",p));
104 17 : PetscFunctionReturn(PETSC_SUCCESS);
105 : }
106 13 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
107 13 : if (ds->compact) {
108 1 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
109 1 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
110 1 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
111 1 : rowsa = n;
112 1 : colsa = ds->extrarow? m+1: m;
113 1 : rowsb = p;
114 1 : colsb = ds->extrarow? m+1: m;
115 1 : if (format == PETSC_VIEWER_ASCII_MATLAB) {
116 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rowsa,colsa));
117 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",2*ds->n));
118 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
119 0 : 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 0 : 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 0 : if (n>m) { /* A lower bidiagonal */
122 0 : 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 0 : 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 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_T]));
127 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rowsb,colsb));
128 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",2*ds->n));
129 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
130 0 : for (i=0;i<rowsb;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)S[i]));
131 0 : for (i=0;i<colsb-1;i++) {
132 0 : r = PetscMax(i+2,ds->k+1);
133 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,r,(double)T[i+2*ds->ld]));
134 : }
135 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_D]));
136 : } else {
137 1 : PetscCall(PetscViewerASCIIPrintf(viewer,"Matrix %s =\n",DSMatName[DS_MAT_T]));
138 11 : for (i=0;i<rowsa;i++) {
139 110 : for (j=0;j<colsa;j++) {
140 100 : if (i==j) value = T[i];
141 90 : else if (i<ds->k && j==ds->k) value = T[i+ds->ld];
142 90 : else if (n>m && i==j+1 && i>ds->k) value = T[j+ds->ld];
143 90 : else if (n<=m && i+1==j && i>=ds->k) value = T[i+ds->ld];
144 : else value = 0.0;
145 100 : PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
146 : }
147 10 : PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
148 : }
149 1 : PetscCall(PetscViewerASCIIPrintf(viewer,"Matrix %s =\n",DSMatName[DS_MAT_D]));
150 11 : for (i=0;i<rowsb;i++) {
151 110 : for (j=0;j<colsb;j++) {
152 100 : if (i==j) value = S[i];
153 90 : else if (i<ds->k && j==ds->k) value = T[PetscMin(i,j)+2*ds->ld];
154 90 : else if (i+1==j && i>=ds->k) value = T[i+2*ds->ld];
155 : else value = 0.0;
156 100 : PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
157 : }
158 10 : PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
159 : }
160 : }
161 1 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
162 1 : PetscCall(PetscViewerFlush(viewer));
163 1 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
164 1 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
165 : } else {
166 12 : PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
167 12 : PetscCall(DSViewMat(ds,viewer,DS_MAT_B));
168 : }
169 13 : if (ds->state>DS_STATE_INTERMEDIATE) {
170 0 : PetscCall(DSViewMat(ds,viewer,DS_MAT_X));
171 0 : PetscCall(DSViewMat(ds,viewer,DS_MAT_U));
172 0 : PetscCall(DSViewMat(ds,viewer,DS_MAT_V));
173 : }
174 13 : PetscFunctionReturn(PETSC_SUCCESS);
175 : }
176 :
177 17 : static PetscErrorCode DSVectors_GSVD(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
178 : {
179 17 : PetscFunctionBegin;
180 17 : switch (mat) {
181 0 : case DS_MAT_U:
182 : case DS_MAT_V:
183 0 : if (rnorm) *rnorm = 0.0;
184 : break;
185 : case DS_MAT_X:
186 : break;
187 0 : default:
188 0 : SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
189 : }
190 17 : PetscFunctionReturn(PETSC_SUCCESS);
191 : }
192 :
193 644 : static PetscErrorCode DSSort_GSVD(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
194 : {
195 644 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
196 644 : PetscInt t,l,ld=ds->ld,i,*perm,*perm2;
197 644 : PetscReal *T=NULL,*D=NULL,*eig;
198 644 : PetscScalar *A=NULL,*B=NULL;
199 644 : PetscBool compact=ds->compact;
200 :
201 644 : PetscFunctionBegin;
202 644 : if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
203 644 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
204 644 : l = ds->l;
205 644 : t = ds->t;
206 644 : perm = ds->perm;
207 644 : PetscCall(PetscMalloc2(t,&eig,t,&perm2));
208 644 : if (compact) {
209 625 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
210 625 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
211 6623 : for (i=0;i<t;i++) eig[i] = (D[i]==0)?PETSC_INFINITY:T[i]/D[i];
212 : } else {
213 19 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
214 19 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
215 230 : 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 644 : PetscCall(DSSortEigenvaluesReal_Private(ds,eig,perm));
218 644 : PetscCall(PetscArraycpy(perm2,perm,t));
219 6612 : for (i=l;i<t;i++) wr[i] = eig[perm[i]];
220 644 : if (compact) {
221 625 : PetscCall(PetscArraycpy(eig,T,t));
222 6382 : for (i=l;i<t;i++) T[i] = eig[perm[i]];
223 625 : PetscCall(PetscArraycpy(eig,D,t));
224 6382 : for (i=l;i<t;i++) D[i] = eig[perm[i]];
225 625 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
226 625 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
227 : } else {
228 230 : for (i=l;i<t;i++) eig[i] = PetscRealPart(A[i*(1+ld)]);
229 230 : for (i=l;i<t;i++) A[i*(1+ld)] = eig[perm[i]];
230 230 : for (i=l;i<t;i++) eig[i] = PetscRealPart(B[i*(1+ld)]);
231 230 : for (i=l;i<t;i++) B[i*(1+ld)] = eig[perm[i]];
232 19 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
233 19 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
234 : }
235 644 : PetscCall(DSPermuteColumns_Private(ds,l,t,ds->n,DS_MAT_U,perm2));
236 644 : PetscCall(PetscArraycpy(perm2,perm,t));
237 644 : PetscCall(DSPermuteColumns_Private(ds,l,t,ctx->m,DS_MAT_X,perm2));
238 644 : PetscCall(DSPermuteColumns_Private(ds,l,t,ctx->p,DS_MAT_V,perm));
239 644 : PetscCall(PetscFree2(eig,perm2));
240 644 : PetscFunctionReturn(PETSC_SUCCESS);
241 : }
242 :
243 622 : static PetscErrorCode DSUpdateExtraRow_GSVD(DS ds)
244 : {
245 622 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
246 622 : PetscInt i;
247 622 : PetscBLASInt n=0,m=0,ld=0;
248 622 : const PetscScalar *U,*V;
249 622 : PetscReal *T,*e,*f,alpha,beta,betah;
250 :
251 622 : PetscFunctionBegin;
252 622 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
253 622 : PetscCheck(ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-compact storage");
254 622 : PetscCall(PetscBLASIntCast(ds->n,&n));
255 622 : PetscCall(PetscBLASIntCast(ctx->m,&m));
256 622 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
257 622 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
258 622 : e = T+ld;
259 622 : f = T+2*ld;
260 622 : PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_U],&U));
261 622 : PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_V],&V));
262 622 : if (n<=m) { /* upper variant, A is square upper bidiagonal */
263 248 : beta = e[m-1]; /* in compact, we assume all entries are zero except the last one */
264 248 : betah = f[m-1];
265 2475 : for (i=0;i<m;i++) {
266 2227 : e[i] = PetscRealPart(beta*U[m-1+i*ld]);
267 2227 : f[i] = PetscRealPart(betah*V[m-1+i*ld]);
268 : }
269 : } else { /* lower variant, A is (m+1)xm lower bidiagonal */
270 374 : alpha = T[m];
271 374 : betah = f[m-1];
272 4115 : for (i=0;i<m;i++) {
273 3741 : e[i] = PetscRealPart(alpha*U[m+i*ld]);
274 3741 : f[i] = PetscRealPart(betah*V[m-1+i*ld]);
275 : }
276 374 : T[m] = PetscRealPart(alpha*U[m+m*ld]);
277 : }
278 622 : ds->k = m;
279 622 : PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_U],&U));
280 622 : PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_V],&V));
281 622 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
282 622 : PetscFunctionReturn(PETSC_SUCCESS);
283 : }
284 :
285 606 : static PetscErrorCode DSTruncate_GSVD(DS ds,PetscInt n,PetscBool trim)
286 : {
287 606 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
288 606 : PetscScalar *U;
289 606 : PetscReal *T;
290 606 : PetscInt i,m=ctx->m,ld=ds->ld;
291 606 : PetscBool lower=(ds->n>ctx->m)?PETSC_TRUE:PETSC_FALSE;
292 :
293 606 : PetscFunctionBegin;
294 606 : PetscCheck(ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-compact storage");
295 606 : if (trim) {
296 54 : ds->l = 0;
297 54 : ds->k = 0;
298 54 : ds->n = lower? n+1: n;
299 54 : ctx->m = n;
300 54 : ctx->p = n;
301 54 : ds->t = ds->n; /* truncated length equal to the new dimension */
302 54 : ctx->tm = ctx->m; /* must also keep the previous dimension of X */
303 54 : ctx->tp = ctx->p; /* must also keep the previous dimension of V */
304 : } else {
305 552 : if (lower) {
306 : /* move value of diagonal element of arrow (alpha) */
307 326 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
308 326 : T[n] = T[m];
309 326 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
310 : /* copy last column of U so that it updates the next initial vector of U1 */
311 326 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
312 3856 : for (i=0;i<=m;i++) U[i+n*ld] = U[i+m*ld];
313 326 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
314 : }
315 552 : ds->k = (ds->extrarow)? n: 0;
316 552 : ds->t = ds->n; /* truncated length equal to previous dimension */
317 552 : ctx->tm = ctx->m; /* must also keep the previous dimension of X */
318 552 : ctx->tp = ctx->p; /* must also keep the previous dimension of V */
319 552 : ds->n = lower? n+1: n;
320 552 : ctx->m = n;
321 552 : ctx->p = n;
322 : }
323 606 : PetscFunctionReturn(PETSC_SUCCESS);
324 : }
325 :
326 637 : static PetscErrorCode DSSwitchFormat_GSVD(DS ds)
327 : {
328 637 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
329 637 : PetscReal *T,*D;
330 637 : PetscScalar *A,*B;
331 637 : PetscInt i,n=ds->n,k=ds->k,ld=ds->ld,m=ctx->m;
332 :
333 637 : PetscFunctionBegin;
334 637 : 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 637 : PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_A],&A));
338 637 : PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_B],&B));
339 637 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
340 637 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
341 637 : PetscCall(PetscArrayzero(A,ld*ld));
342 637 : PetscCall(PetscArrayzero(B,ld*ld));
343 3120 : for (i=0;i<k;i++) {
344 2483 : A[i+i*ld] = T[i];
345 2483 : A[i+k*ld] = T[i+ld];
346 2483 : B[i+i*ld] = D[i];
347 2483 : B[i+k*ld] = T[i+2*ld];
348 : }
349 : /* B is upper bidiagonal */
350 637 : B[k+k*ld] = D[k];
351 3612 : for (i=k+1;i<m;i++) {
352 2975 : B[i+i*ld] = D[i];
353 2975 : B[i-1+i*ld] = T[i-1+2*ld];
354 : }
355 : /* A can be upper (square) or lower bidiagonal */
356 4242 : for (i=k;i<m;i++) A[i+i*ld] = T[i];
357 2873 : if (n>m) for (i=k;i<m;i++) A[i+1+i*ld] = T[i+ld];
358 1369 : else for (i=k+1;i<m;i++) A[i-1+i*ld] = T[i-1+ld];
359 637 : PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_A],&A));
360 637 : PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_B],&B));
361 637 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
362 637 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
363 637 : PetscFunctionReturn(PETSC_SUCCESS);
364 : }
365 :
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 644 : static PetscErrorCode DSSolve_GSVD(DS ds,PetscScalar *wr,PetscScalar *wi)
371 : {
372 644 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
373 644 : PetscInt i,j;
374 644 : PetscBLASInt n1,m1,info,lc = 0,n = 0,m = 0,p = 0,p1,l,k,q,ld,off,lwork,r;
375 644 : PetscScalar *A,*B,*X,*U,*V,sone=1.0,smone=-1.0;
376 644 : PetscReal *alpha,*beta,*T,*D;
377 : #if !defined(SLEPC_MISSING_LAPACK_GGSVD3)
378 644 : PetscScalar a,dummy;
379 644 : PetscReal rdummy;
380 644 : PetscBLASInt idummy;
381 : #endif
382 :
383 644 : PetscFunctionBegin;
384 644 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
385 644 : PetscCall(PetscBLASIntCast(ds->n,&m));
386 644 : PetscCall(PetscBLASIntCast(ctx->m,&n));
387 644 : PetscCall(PetscBLASIntCast(ctx->p,&p));
388 644 : PetscCall(PetscBLASIntCast(ds->l,&lc));
389 644 : 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 644 : PetscCheck(!ds->compact || (p==n && (m==p || m==p+1)),PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Dimensions not supported in compact format");
392 644 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
393 644 : n1 = n-lc; /* n1 = size of leading block, excl. locked + size of trailing block */
394 644 : m1 = m-lc;
395 644 : p1 = p-lc;
396 644 : off = lc+lc*ld;
397 644 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
398 644 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
399 644 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
400 644 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
401 644 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
402 644 : PetscCall(PetscArrayzero(X,ld*ld));
403 885 : for (i=0;i<lc;i++) X[i+i*ld] = 1.0;
404 644 : PetscCall(PetscArrayzero(U,ld*ld));
405 885 : for (i=0;i<lc;i++) U[i+i*ld] = 1.0;
406 644 : PetscCall(PetscArrayzero(V,ld*ld));
407 885 : for (i=0;i<lc;i++) V[i+i*ld] = 1.0;
408 644 : if (ds->compact) PetscCall(DSSwitchFormat_GSVD(ds));
409 :
410 : #if !defined(SLEPC_MISSING_LAPACK_GGSVD3)
411 : /* workspace query and memory allocation */
412 644 : lwork = -1;
413 : #if !defined (PETSC_USE_COMPLEX)
414 644 : 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 644 : 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
420 :
421 : #if !defined (PETSC_USE_COMPLEX)
422 644 : PetscCall(DSAllocateWork_Private(ds,lwork,2*ds->ld,ds->ld));
423 644 : alpha = ds->rwork;
424 644 : beta = ds->rwork+ds->ld;
425 644 : 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 644 : SlepcCheckLapackInfo("ggsvd3",info);
433 :
434 : #else /* defined(SLEPC_MISSING_LAPACK_GGSVD3) */
435 :
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);
449 :
450 : #endif
451 :
452 644 : PetscCheck(k+l>=n1,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"The rank deficient case not supported yet");
453 644 : if (ds->compact) {
454 625 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
455 625 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
456 : /* R is the identity matrix (except the sign) */
457 6382 : for (i=lc;i<n;i++) {
458 5757 : if (PetscRealPart(A[i+i*ld])<0.0) { /* scale column i */
459 37340 : for (j=lc;j<n;j++) X[j+i*ld] = -X[j+i*ld];
460 : }
461 : }
462 625 : PetscCall(PetscArrayzero(T+ld,m-1));
463 625 : PetscCall(PetscArrayzero(T+2*ld,n-1));
464 6382 : for (i=lc;i<n;i++) {
465 5757 : T[i] = alpha[i-lc];
466 5757 : D[i] = beta[i-lc];
467 5757 : if (D[i]==0.0) wr[i] = PETSC_INFINITY;
468 5757 : else wr[i] = T[i]/D[i];
469 : }
470 625 : ds->t = n;
471 625 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
472 625 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
473 : } else {
474 : /* X = X*inv(R) */
475 19 : q = PetscMin(m,n);
476 19 : PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n,&q,&sone,A,&ld,X,&ld));
477 19 : if (m<n) {
478 4 : r = n-m;
479 4 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&r,&m,&sone,X,&ld,A,&ld,&smone,X+m*ld,&ld));
480 4 : PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n,&r,&sone,B+m*ld,&ld,X+m*ld,&ld));
481 : }
482 19 : if (k>0) {
483 13 : for (i=k;i<PetscMin(m,k+l);i++) {
484 11 : PetscCall(PetscArraycpy(X+(i-k)*ld,X+i*ld,ld));
485 11 : PetscCall(PetscArraycpy(U+(i-k)*ld,U+i*ld,ld));
486 : }
487 : }
488 : /* singular values */
489 19 : PetscCall(PetscArrayzero(A,ld*ld));
490 19 : PetscCall(PetscArrayzero(B,ld*ld));
491 230 : for (j=k;j<PetscMin(m,k+l);j++) {
492 211 : A[(j-k)*(1+ld)] = alpha[j];
493 211 : B[(j-k)*(1+ld)] = beta[j];
494 211 : wr[j-k] = alpha[j]/beta[j];
495 : }
496 19 : ds->t = PetscMin(m,k+l)-k; /* set number of computed values */
497 : }
498 644 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
499 644 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
500 644 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
501 644 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
502 644 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
503 644 : PetscFunctionReturn(PETSC_SUCCESS);
504 : }
505 :
506 24 : static PetscErrorCode DSCond_GSVD(DS ds,PetscReal *cond)
507 : {
508 24 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
509 24 : PetscBLASInt lwork,lrwork=0,info,m,n,p,ld;
510 24 : PetscScalar *A,*work;
511 24 : const PetscScalar *M;
512 24 : PetscReal *sigma,conda,condb;
513 : #if defined(PETSC_USE_COMPLEX)
514 : PetscReal *rwork;
515 : #endif
516 :
517 24 : PetscFunctionBegin;
518 24 : PetscCall(PetscBLASIntCast(ds->n,&m));
519 24 : PetscCall(PetscBLASIntCast(ctx->m,&n));
520 24 : PetscCall(PetscBLASIntCast(ctx->p,&p));
521 24 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
522 24 : lwork = 5*n;
523 : #if defined(PETSC_USE_COMPLEX)
524 : lrwork = 5*n;
525 : #endif
526 24 : PetscCall(DSAllocateWork_Private(ds,ld*n+lwork,n+lrwork,0));
527 24 : A = ds->work;
528 24 : work = ds->work+ld*n;
529 24 : sigma = ds->rwork;
530 : #if defined(PETSC_USE_COMPLEX)
531 : rwork = ds->rwork+n;
532 : #endif
533 24 : if (ds->compact) PetscCall(DSSwitchFormat_GSVD(ds));
534 :
535 24 : PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&M));
536 24 : PetscCall(PetscArraycpy(A,M,ld*n));
537 24 : 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 24 : PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,&info));
542 : #endif
543 24 : SlepcCheckLapackInfo("gesvd",info);
544 24 : conda = sigma[0]/sigma[PetscMin(m,n)-1];
545 :
546 24 : PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&M));
547 24 : PetscCall(PetscArraycpy(A,M,ld*n));
548 24 : 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 24 : PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&p,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,&info));
553 : #endif
554 24 : SlepcCheckLapackInfo("gesvd",info);
555 24 : condb = sigma[0]/sigma[PetscMin(p,n)-1];
556 :
557 24 : *cond = PetscMax(conda,condb);
558 24 : PetscFunctionReturn(PETSC_SUCCESS);
559 : }
560 :
561 : #if !defined(PETSC_HAVE_MPIUNI)
562 29 : static PetscErrorCode DSSynchronize_GSVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
563 : {
564 29 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
565 29 : PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
566 29 : PetscMPIInt m=ctx->m,rank,off=0,size,n,ldn,ld3;
567 29 : PetscScalar *A,*U,*V,*X;
568 29 : PetscReal *T;
569 :
570 29 : PetscFunctionBegin;
571 29 : if (ds->compact) kr = 3*ld;
572 5 : else k = 2*(m-l)*ld;
573 29 : if (ds->state>DS_STATE_RAW) k += 3*(m-l)*ld;
574 29 : if (eigr) k += m-l;
575 29 : PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
576 29 : PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
577 29 : PetscCall(PetscMPIIntCast(m-l,&n));
578 29 : PetscCall(PetscMPIIntCast(ld*(m-l),&ldn));
579 29 : PetscCall(PetscMPIIntCast(3*ld,&ld3));
580 29 : if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
581 5 : else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
582 29 : if (ds->state>DS_STATE_RAW) {
583 29 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
584 29 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
585 29 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
586 : }
587 29 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
588 29 : if (!rank) {
589 14 : if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
590 2 : else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
591 14 : if (ds->state>DS_STATE_RAW) {
592 14 : PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
593 14 : PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
594 14 : PetscCallMPI(MPI_Pack(X+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
595 : }
596 14 : if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
597 : }
598 58 : PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
599 29 : if (rank) {
600 15 : if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
601 3 : else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
602 15 : if (ds->state>DS_STATE_RAW) {
603 15 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
604 15 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
605 15 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,X+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
606 : }
607 15 : if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
608 : }
609 29 : if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
610 5 : else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
611 29 : if (ds->state>DS_STATE_RAW) {
612 29 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
613 29 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
614 29 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
615 : }
616 29 : PetscFunctionReturn(PETSC_SUCCESS);
617 : }
618 : #endif
619 :
620 1907 : static PetscErrorCode DSMatGetSize_GSVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
621 : {
622 1907 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
623 :
624 1907 : PetscFunctionBegin;
625 1907 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
626 1907 : switch (t) {
627 19 : case DS_MAT_A:
628 19 : *rows = ds->n;
629 19 : *cols = ds->extrarow? ctx->m+1: ctx->m;
630 19 : break;
631 19 : case DS_MAT_B:
632 19 : *rows = ctx->p;
633 19 : *cols = ds->extrarow? ctx->m+1: ctx->m;
634 19 : break;
635 0 : case DS_MAT_T:
636 0 : *rows = ds->n;
637 0 : *cols = PetscDefined(USE_COMPLEX)? 2: 3;
638 0 : break;
639 0 : case DS_MAT_D:
640 0 : *rows = ctx->p;
641 0 : *cols = 1;
642 0 : break;
643 623 : case DS_MAT_U:
644 623 : *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
645 623 : *cols = ds->n;
646 623 : break;
647 623 : case DS_MAT_V:
648 623 : *rows = ds->state==DS_STATE_TRUNCATED? ctx->tp: ctx->p;
649 623 : *cols = ctx->p;
650 623 : break;
651 623 : case DS_MAT_X:
652 623 : *rows = ds->state==DS_STATE_TRUNCATED? ctx->tm: ctx->m;
653 623 : *cols = ctx->m;
654 623 : break;
655 0 : default:
656 0 : SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
657 : }
658 1907 : PetscFunctionReturn(PETSC_SUCCESS);
659 : }
660 :
661 644 : static PetscErrorCode DSGSVDSetDimensions_GSVD(DS ds,PetscInt m,PetscInt p)
662 : {
663 644 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
664 :
665 644 : PetscFunctionBegin;
666 644 : DSCheckAlloc(ds,1);
667 644 : if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
668 0 : ctx->m = ds->ld;
669 : } else {
670 644 : PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
671 644 : ctx->m = m;
672 : }
673 644 : if (p==PETSC_DECIDE || p==PETSC_DEFAULT) {
674 3 : ctx->p = ds->n;
675 : } else {
676 641 : PetscCheck(p>0 && p<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of p. Must be between 1 and ld");
677 641 : ctx->p = p;
678 : }
679 644 : PetscFunctionReturn(PETSC_SUCCESS);
680 : }
681 :
682 : /*@
683 : DSGSVDSetDimensions - Sets the number of columns and rows for a DSGSVD.
684 :
685 : Logically Collective
686 :
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)
691 :
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().
696 :
697 : Level: intermediate
698 :
699 : .seealso: DSGSVDGetDimensions(), DSSetDimensions()
700 : @*/
701 644 : PetscErrorCode DSGSVDSetDimensions(DS ds,PetscInt m,PetscInt p)
702 : {
703 644 : PetscFunctionBegin;
704 644 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
705 2576 : PetscValidLogicalCollectiveInt(ds,m,2);
706 2576 : PetscValidLogicalCollectiveInt(ds,p,3);
707 644 : PetscTryMethod(ds,"DSGSVDSetDimensions_C",(DS,PetscInt,PetscInt),(ds,m,p));
708 644 : PetscFunctionReturn(PETSC_SUCCESS);
709 : }
710 :
711 12 : static PetscErrorCode DSGSVDGetDimensions_GSVD(DS ds,PetscInt *m,PetscInt *p)
712 : {
713 12 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
714 :
715 12 : PetscFunctionBegin;
716 12 : if (m) *m = ctx->m;
717 12 : if (p) *p = ctx->p;
718 12 : PetscFunctionReturn(PETSC_SUCCESS);
719 : }
720 :
721 : /*@
722 : DSGSVDGetDimensions - Returns the number of columns and rows for a DSGSVD.
723 :
724 : Not Collective
725 :
726 : Input Parameter:
727 : . ds - the direct solver context
728 :
729 : Output Parameters:
730 : + m - the number of columns
731 : - p - the number of rows for the second matrix (B)
732 :
733 : Level: intermediate
734 :
735 : .seealso: DSGSVDSetDimensions()
736 : @*/
737 12 : PetscErrorCode DSGSVDGetDimensions(DS ds,PetscInt *m,PetscInt *p)
738 : {
739 12 : PetscFunctionBegin;
740 12 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
741 12 : PetscUseMethod(ds,"DSGSVDGetDimensions_C",(DS,PetscInt*,PetscInt*),(ds,m,p));
742 12 : PetscFunctionReturn(PETSC_SUCCESS);
743 : }
744 :
745 72 : static PetscErrorCode DSDestroy_GSVD(DS ds)
746 : {
747 72 : PetscFunctionBegin;
748 72 : PetscCall(PetscFree(ds->data));
749 72 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDSetDimensions_C",NULL));
750 72 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDGetDimensions_C",NULL));
751 72 : PetscFunctionReturn(PETSC_SUCCESS);
752 : }
753 :
754 : /*MC
755 : DSGSVD - Dense Generalized Singular Value Decomposition.
756 :
757 : Level: beginner
758 :
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.
767 :
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().
771 :
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().
774 :
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.
783 :
784 : Allowed arguments of DSVectors() are DS_MAT_U, DS_MAT_V and DS_MAT_X.
785 :
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)
791 :
792 : Implemented methods:
793 : . 0 - Lapack (_ggsvd3 if available, or _ggsvd)
794 :
795 : .seealso: DSCreate(), DSSetType(), DSType, DSGSVDSetDimensions()
796 : M*/
797 72 : SLEPC_EXTERN PetscErrorCode DSCreate_GSVD(DS ds)
798 : {
799 72 : DS_GSVD *ctx;
800 :
801 72 : PetscFunctionBegin;
802 72 : PetscCall(PetscNew(&ctx));
803 72 : ds->data = (void*)ctx;
804 :
805 72 : ds->ops->allocate = DSAllocate_GSVD;
806 72 : ds->ops->view = DSView_GSVD;
807 72 : ds->ops->vectors = DSVectors_GSVD;
808 72 : ds->ops->sort = DSSort_GSVD;
809 72 : ds->ops->solve[0] = DSSolve_GSVD;
810 : #if !defined(PETSC_HAVE_MPIUNI)
811 72 : ds->ops->synchronize = DSSynchronize_GSVD;
812 : #endif
813 72 : ds->ops->truncate = DSTruncate_GSVD;
814 72 : ds->ops->update = DSUpdateExtraRow_GSVD;
815 72 : ds->ops->cond = DSCond_GSVD;
816 72 : ds->ops->matgetsize = DSMatGetSize_GSVD;
817 72 : ds->ops->destroy = DSDestroy_GSVD;
818 72 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDSetDimensions_C",DSGSVDSetDimensions_GSVD));
819 72 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDGetDimensions_C",DSGSVDGetDimensions_GSVD));
820 72 : PetscFunctionReturn(PETSC_SUCCESS);
821 : }
|
