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 66 : static PetscErrorCode DSAllocate_GSVD(DS ds,PetscInt ld)
21 : {
22 66 : PetscFunctionBegin;
23 66 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
24 66 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
25 66 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
26 66 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
27 66 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_V));
28 66 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
29 66 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_D));
30 66 : PetscCall(PetscFree(ds->perm));
31 66 : PetscCall(PetscMalloc1(ld,&ds->perm));
32 66 : 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 560 : static PetscErrorCode DSSort_GSVD(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
194 : {
195 560 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
196 560 : PetscInt t,l,ld=ds->ld,i,*perm,*perm2;
197 560 : PetscReal *T=NULL,*D=NULL,*eig;
198 560 : PetscScalar *A=NULL,*B=NULL;
199 560 : PetscBool compact=ds->compact;
200 :
201 560 : PetscFunctionBegin;
202 560 : if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
203 560 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
204 560 : l = ds->l;
205 560 : t = ds->t;
206 560 : perm = ds->perm;
207 560 : PetscCall(PetscMalloc2(t,&eig,t,&perm2));
208 560 : if (compact) {
209 541 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
210 541 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
211 5465 : 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 560 : PetscCall(DSSortEigenvaluesReal_Private(ds,eig,perm));
218 560 : PetscCall(PetscArraycpy(perm2,perm,t));
219 5492 : for (i=l;i<t;i++) wr[i] = eig[perm[i]];
220 560 : if (compact) {
221 541 : PetscCall(PetscArraycpy(eig,T,t));
222 5262 : for (i=l;i<t;i++) T[i] = eig[perm[i]];
223 541 : PetscCall(PetscArraycpy(eig,D,t));
224 5262 : for (i=l;i<t;i++) D[i] = eig[perm[i]];
225 541 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
226 541 : 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 560 : PetscCall(DSPermuteColumns_Private(ds,l,t,ds->n,DS_MAT_U,perm2));
236 560 : PetscCall(PetscArraycpy(perm2,perm,t));
237 560 : PetscCall(DSPermuteColumns_Private(ds,l,t,ctx->m,DS_MAT_X,perm2));
238 560 : PetscCall(DSPermuteColumns_Private(ds,l,t,ctx->p,DS_MAT_V,perm));
239 560 : PetscCall(PetscFree2(eig,perm2));
240 560 : PetscFunctionReturn(PETSC_SUCCESS);
241 : }
242 :
243 538 : static PetscErrorCode DSUpdateExtraRow_GSVD(DS ds)
244 : {
245 538 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
246 538 : PetscInt i;
247 538 : PetscBLASInt n=0,m=0,ld=0;
248 538 : const PetscScalar *U,*V;
249 538 : PetscReal *T,*e,*f,alpha,beta,betah;
250 :
251 538 : PetscFunctionBegin;
252 538 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
253 538 : PetscCheck(ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-compact storage");
254 538 : PetscCall(PetscBLASIntCast(ds->n,&n));
255 538 : PetscCall(PetscBLASIntCast(ctx->m,&m));
256 538 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
257 538 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
258 538 : e = T+ld;
259 538 : f = T+2*ld;
260 538 : PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_U],&U));
261 538 : PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_V],&V));
262 538 : if (n<=m) { /* upper variant, A is square upper bidiagonal */
263 224 : beta = e[m-1]; /* in compact, we assume all entries are zero except the last one */
264 224 : betah = f[m-1];
265 2209 : for (i=0;i<m;i++) {
266 1985 : e[i] = PetscRealPart(beta*U[m-1+i*ld]);
267 1985 : f[i] = PetscRealPart(betah*V[m-1+i*ld]);
268 : }
269 : } else { /* lower variant, A is (m+1)xm lower bidiagonal */
270 314 : alpha = T[m];
271 314 : betah = f[m-1];
272 3223 : for (i=0;i<m;i++) {
273 2909 : e[i] = PetscRealPart(alpha*U[m+i*ld]);
274 2909 : f[i] = PetscRealPart(betah*V[m-1+i*ld]);
275 : }
276 314 : T[m] = PetscRealPart(alpha*U[m+m*ld]);
277 : }
278 538 : ds->k = m;
279 538 : PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_U],&U));
280 538 : PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_V],&V));
281 538 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
282 538 : PetscFunctionReturn(PETSC_SUCCESS);
283 : }
284 :
285 529 : static PetscErrorCode DSTruncate_GSVD(DS ds,PetscInt n,PetscBool trim)
286 : {
287 529 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
288 529 : PetscScalar *U;
289 529 : PetscReal *T;
290 529 : PetscInt i,m=ctx->m,ld=ds->ld;
291 529 : PetscBool lower=(ds->n>ctx->m)?PETSC_TRUE:PETSC_FALSE;
292 :
293 529 : PetscFunctionBegin;
294 529 : PetscCheck(ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-compact storage");
295 529 : if (trim) {
296 48 : ds->l = 0;
297 48 : ds->k = 0;
298 48 : ds->n = lower? n+1: n;
299 48 : ctx->m = n;
300 48 : ctx->p = n;
301 48 : ds->t = ds->n; /* truncated length equal to the new dimension */
302 48 : ctx->tm = ctx->m; /* must also keep the previous dimension of X */
303 48 : ctx->tp = ctx->p; /* must also keep the previous dimension of V */
304 : } else {
305 481 : if (lower) {
306 : /* move value of diagonal element of arrow (alpha) */
307 275 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
308 275 : T[n] = T[m];
309 275 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
310 : /* copy last column of U so that it updates the next initial vector of U1 */
311 275 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
312 3084 : for (i=0;i<=m;i++) U[i+n*ld] = U[i+m*ld];
313 275 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
314 : }
315 481 : ds->k = (ds->extrarow)? n: 0;
316 481 : ds->t = ds->n; /* truncated length equal to previous dimension */
317 481 : ctx->tm = ctx->m; /* must also keep the previous dimension of X */
318 481 : ctx->tp = ctx->p; /* must also keep the previous dimension of V */
319 481 : ds->n = lower? n+1: n;
320 481 : ctx->m = n;
321 481 : ctx->p = n;
322 : }
323 529 : PetscFunctionReturn(PETSC_SUCCESS);
324 : }
325 :
326 553 : static PetscErrorCode DSSwitchFormat_GSVD(DS ds)
327 : {
328 553 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
329 553 : PetscReal *T,*D;
330 553 : PetscScalar *A,*B;
331 553 : PetscInt i,n=ds->n,k=ds->k,ld=ds->ld,m=ctx->m;
332 :
333 553 : PetscFunctionBegin;
334 553 : 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 553 : PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_A],&A));
338 553 : PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_B],&B));
339 553 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
340 553 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
341 553 : PetscCall(PetscArrayzero(A,ld*ld));
342 553 : PetscCall(PetscArrayzero(B,ld*ld));
343 2587 : for (i=0;i<k;i++) {
344 2034 : A[i+i*ld] = T[i];
345 2034 : A[i+k*ld] = T[i+ld];
346 2034 : B[i+i*ld] = D[i];
347 2034 : B[i+k*ld] = T[i+2*ld];
348 : }
349 : /* B is upper bidiagonal */
350 553 : B[k+k*ld] = D[k];
351 2987 : for (i=k+1;i<m;i++) {
352 2434 : B[i+i*ld] = D[i];
353 2434 : B[i-1+i*ld] = T[i-1+2*ld];
354 : }
355 : /* A can be upper (square) or lower bidiagonal */
356 3533 : for (i=k;i<m;i++) A[i+i*ld] = T[i];
357 2331 : if (n>m) for (i=k;i<m;i++) A[i+1+i*ld] = T[i+ld];
358 1202 : else for (i=k+1;i<m;i++) A[i-1+i*ld] = T[i-1+ld];
359 553 : PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_A],&A));
360 553 : PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_B],&B));
361 553 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
362 553 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
363 553 : 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 560 : static PetscErrorCode DSSolve_GSVD(DS ds,PetscScalar *wr,PetscScalar *wi)
371 : {
372 560 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
373 560 : PetscInt i,j;
374 560 : PetscBLASInt n1,m1,info,lc = 0,n = 0,m = 0,p = 0,p1,l,k,q,ld,off,lwork,r;
375 560 : PetscScalar *A,*B,*X,*U,*V,sone=1.0,smone=-1.0;
376 560 : PetscReal *alpha,*beta,*T,*D;
377 : #if !defined(SLEPC_MISSING_LAPACK_GGSVD3)
378 560 : PetscScalar a,dummy;
379 560 : PetscReal rdummy;
380 560 : PetscBLASInt idummy;
381 : #endif
382 :
383 560 : PetscFunctionBegin;
384 560 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
385 560 : PetscCall(PetscBLASIntCast(ds->n,&m));
386 560 : PetscCall(PetscBLASIntCast(ctx->m,&n));
387 560 : PetscCall(PetscBLASIntCast(ctx->p,&p));
388 560 : PetscCall(PetscBLASIntCast(ds->l,&lc));
389 560 : 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 560 : PetscCheck(!ds->compact || (p==n && (m==p || m==p+1)),PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Dimensions not supported in compact format");
392 560 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
393 560 : n1 = n-lc; /* n1 = size of leading block, excl. locked + size of trailing block */
394 560 : m1 = m-lc;
395 560 : p1 = p-lc;
396 560 : off = lc+lc*ld;
397 560 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
398 560 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
399 560 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
400 560 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
401 560 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
402 560 : PetscCall(PetscArrayzero(X,ld*ld));
403 763 : for (i=0;i<lc;i++) X[i+i*ld] = 1.0;
404 560 : PetscCall(PetscArrayzero(U,ld*ld));
405 763 : for (i=0;i<lc;i++) U[i+i*ld] = 1.0;
406 560 : PetscCall(PetscArrayzero(V,ld*ld));
407 763 : for (i=0;i<lc;i++) V[i+i*ld] = 1.0;
408 560 : if (ds->compact) PetscCall(DSSwitchFormat_GSVD(ds));
409 :
410 : #if !defined(SLEPC_MISSING_LAPACK_GGSVD3)
411 : /* workspace query and memory allocation */
412 560 : 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 560 : 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 560 : PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
419 : #endif
420 :
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 560 : PetscCall(DSAllocateWork_Private(ds,lwork,4*ds->ld,ds->ld));
428 560 : alpha = ds->rwork+2*ds->ld;
429 560 : beta = ds->rwork+3*ds->ld;
430 560 : 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 560 : 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 560 : PetscCheck(k+l>=n1,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"The rank deficient case not supported yet");
453 560 : if (ds->compact) {
454 541 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
455 541 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
456 : /* R is the identity matrix (except the sign) */
457 5262 : for (i=lc;i<n;i++) {
458 4721 : if (PetscRealPart(A[i+i*ld])<0.0) { /* scale column i */
459 28742 : for (j=lc;j<n;j++) X[j+i*ld] = -X[j+i*ld];
460 : }
461 : }
462 541 : PetscCall(PetscArrayzero(T+ld,m-1));
463 541 : PetscCall(PetscArrayzero(T+2*ld,n-1));
464 5262 : for (i=lc;i<n;i++) {
465 4721 : T[i] = alpha[i-lc];
466 4721 : D[i] = beta[i-lc];
467 4721 : if (D[i]==0.0) wr[i] = PETSC_INFINITY;
468 4721 : else wr[i] = T[i]/D[i];
469 : }
470 541 : ds->t = n;
471 541 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
472 541 : 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 560 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
499 560 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
500 560 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
501 560 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
502 560 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
503 560 : 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 24 : 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 24 : 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 24 : 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 24 : 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 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 24 : 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 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 31 : static PetscErrorCode DSSynchronize_GSVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
563 : {
564 31 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
565 31 : PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
566 31 : PetscMPIInt m,rank,off=0,size,n,ldn,ld3;
567 31 : PetscScalar *A,*U,*V,*X;
568 31 : PetscReal *T;
569 :
570 31 : PetscFunctionBegin;
571 31 : PetscCall(PetscMPIIntCast(ctx->m,&m));
572 31 : if (ds->compact) kr = 3*ld;
573 5 : else k = 2*(m-l)*ld;
574 31 : if (ds->state>DS_STATE_RAW) k += 3*(m-l)*ld;
575 31 : if (eigr) k += m-l;
576 31 : PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
577 31 : PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
578 31 : PetscCall(PetscMPIIntCast(m-l,&n));
579 31 : PetscCall(PetscMPIIntCast(ld*(m-l),&ldn));
580 31 : PetscCall(PetscMPIIntCast(3*ld,&ld3));
581 31 : if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
582 5 : else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
583 31 : if (ds->state>DS_STATE_RAW) {
584 31 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
585 31 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
586 31 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
587 : }
588 31 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
589 31 : if (!rank) {
590 15 : if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
591 2 : else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
592 15 : if (ds->state>DS_STATE_RAW) {
593 15 : PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
594 15 : PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
595 15 : PetscCallMPI(MPI_Pack(X+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
596 : }
597 15 : if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
598 : }
599 62 : PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
600 31 : if (rank) {
601 16 : if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
602 3 : else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
603 16 : if (ds->state>DS_STATE_RAW) {
604 16 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
605 16 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
606 16 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,X+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
607 : }
608 16 : if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
609 : }
610 31 : if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
611 5 : else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
612 31 : if (ds->state>DS_STATE_RAW) {
613 31 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
614 31 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
615 31 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
616 : }
617 31 : PetscFunctionReturn(PETSC_SUCCESS);
618 : }
619 : #endif
620 :
621 1676 : static PetscErrorCode DSMatGetSize_GSVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
622 : {
623 1676 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
624 :
625 1676 : PetscFunctionBegin;
626 1676 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
627 1676 : switch (t) {
628 19 : case DS_MAT_A:
629 19 : *rows = ds->n;
630 19 : *cols = ds->extrarow? ctx->m+1: ctx->m;
631 19 : break;
632 19 : case DS_MAT_B:
633 19 : *rows = ctx->p;
634 19 : *cols = ds->extrarow? ctx->m+1: ctx->m;
635 19 : break;
636 0 : case DS_MAT_T:
637 0 : *rows = ds->n;
638 0 : *cols = PetscDefined(USE_COMPLEX)? 2: 3;
639 0 : break;
640 0 : case DS_MAT_D:
641 0 : *rows = ctx->p;
642 0 : *cols = 1;
643 0 : break;
644 546 : case DS_MAT_U:
645 546 : *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
646 546 : *cols = ds->n;
647 546 : break;
648 546 : case DS_MAT_V:
649 546 : *rows = ds->state==DS_STATE_TRUNCATED? ctx->tp: ctx->p;
650 546 : *cols = ctx->p;
651 546 : break;
652 546 : case DS_MAT_X:
653 546 : *rows = ds->state==DS_STATE_TRUNCATED? ctx->tm: ctx->m;
654 546 : *cols = ctx->m;
655 546 : break;
656 0 : default:
657 0 : SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
658 : }
659 1676 : PetscFunctionReturn(PETSC_SUCCESS);
660 : }
661 :
662 560 : static PetscErrorCode DSGSVDSetDimensions_GSVD(DS ds,PetscInt m,PetscInt p)
663 : {
664 560 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
665 :
666 560 : PetscFunctionBegin;
667 560 : DSCheckAlloc(ds,1);
668 560 : if (m == PETSC_DETERMINE) {
669 0 : ctx->m = ds->ld;
670 560 : } else if (m != PETSC_CURRENT) {
671 560 : PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
672 560 : ctx->m = m;
673 : }
674 560 : if (p == PETSC_DETERMINE) {
675 3 : ctx->p = ds->n;
676 557 : } else if (p != PETSC_CURRENT) {
677 557 : PetscCheck(p>0 && p<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of p. Must be between 1 and ld");
678 557 : ctx->p = p;
679 : }
680 560 : PetscFunctionReturn(PETSC_SUCCESS);
681 : }
682 :
683 : /*@
684 : DSGSVDSetDimensions - Sets the number of columns and rows for a DSGSVD.
685 :
686 : Logically Collective
687 :
688 : Input Parameters:
689 : + ds - the direct solver context
690 : . m - the number of columns
691 : - p - the number of rows for the second matrix (B)
692 :
693 : Notes:
694 : This call is complementary to DSSetDimensions(), to provide two dimensions
695 : that are specific to this DS type. The number of rows for the first matrix (A)
696 : is set by DSSetDimensions().
697 :
698 : Use PETSC_CURRENT to leave any of the values unchanged. Use PETSC_DETERMINE
699 : to set m to the leading dimension and p to the number of columns of B.
700 :
701 : Level: intermediate
702 :
703 : .seealso: DSGSVDGetDimensions(), DSSetDimensions()
704 : @*/
705 560 : PetscErrorCode DSGSVDSetDimensions(DS ds,PetscInt m,PetscInt p)
706 : {
707 560 : PetscFunctionBegin;
708 560 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
709 1680 : PetscValidLogicalCollectiveInt(ds,m,2);
710 1680 : PetscValidLogicalCollectiveInt(ds,p,3);
711 560 : PetscTryMethod(ds,"DSGSVDSetDimensions_C",(DS,PetscInt,PetscInt),(ds,m,p));
712 560 : PetscFunctionReturn(PETSC_SUCCESS);
713 : }
714 :
715 12 : static PetscErrorCode DSGSVDGetDimensions_GSVD(DS ds,PetscInt *m,PetscInt *p)
716 : {
717 12 : DS_GSVD *ctx = (DS_GSVD*)ds->data;
718 :
719 12 : PetscFunctionBegin;
720 12 : if (m) *m = ctx->m;
721 12 : if (p) *p = ctx->p;
722 12 : PetscFunctionReturn(PETSC_SUCCESS);
723 : }
724 :
725 : /*@
726 : DSGSVDGetDimensions - Returns the number of columns and rows for a DSGSVD.
727 :
728 : Not Collective
729 :
730 : Input Parameter:
731 : . ds - the direct solver context
732 :
733 : Output Parameters:
734 : + m - the number of columns
735 : - p - the number of rows for the second matrix (B)
736 :
737 : Level: intermediate
738 :
739 : .seealso: DSGSVDSetDimensions()
740 : @*/
741 12 : PetscErrorCode DSGSVDGetDimensions(DS ds,PetscInt *m,PetscInt *p)
742 : {
743 12 : PetscFunctionBegin;
744 12 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
745 12 : PetscUseMethod(ds,"DSGSVDGetDimensions_C",(DS,PetscInt*,PetscInt*),(ds,m,p));
746 12 : PetscFunctionReturn(PETSC_SUCCESS);
747 : }
748 :
749 66 : static PetscErrorCode DSDestroy_GSVD(DS ds)
750 : {
751 66 : PetscFunctionBegin;
752 66 : PetscCall(PetscFree(ds->data));
753 66 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDSetDimensions_C",NULL));
754 66 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDGetDimensions_C",NULL));
755 66 : PetscFunctionReturn(PETSC_SUCCESS);
756 : }
757 :
758 0 : static PetscErrorCode DSReallocate_GSVD(DS ds,PetscInt ld)
759 : {
760 0 : PetscInt i,*perm=ds->perm;
761 :
762 0 : PetscFunctionBegin;
763 0 : for (i=0;i<DS_NUM_MAT;i++) {
764 0 : if (i!=DS_MAT_A && i!=DS_MAT_B && i!=DS_MAT_X && i!=DS_MAT_U && i!=DS_MAT_V && i!=DS_MAT_T && i!=DS_MAT_D) PetscCall(MatDestroy(&ds->omat[i]));
765 : }
766 :
767 0 : PetscCall(DSReallocateMat_Private(ds,DS_MAT_A,ld));
768 0 : PetscCall(DSReallocateMat_Private(ds,DS_MAT_B,ld));
769 0 : PetscCall(DSReallocateMat_Private(ds,DS_MAT_X,ld));
770 0 : PetscCall(DSReallocateMat_Private(ds,DS_MAT_U,ld));
771 0 : PetscCall(DSReallocateMat_Private(ds,DS_MAT_V,ld));
772 0 : PetscCall(DSReallocateMat_Private(ds,DS_MAT_T,ld));
773 0 : PetscCall(DSReallocateMat_Private(ds,DS_MAT_D,ld));
774 :
775 0 : PetscCall(PetscMalloc1(ld,&ds->perm));
776 0 : PetscCall(PetscArraycpy(ds->perm,perm,ds->ld));
777 0 : PetscCall(PetscFree(perm));
778 0 : PetscFunctionReturn(PETSC_SUCCESS);
779 : }
780 :
781 : /*MC
782 : DSGSVD - Dense Generalized Singular Value Decomposition.
783 :
784 : Level: beginner
785 :
786 : Notes:
787 : The problem is expressed as A*X = U*C, B*X = V*S, where A and B are
788 : matrices with the same number of columns, m, U and V are orthogonal
789 : (unitary), and X is an mxm invertible matrix. The DS object does not
790 : expose matrices C and S, instead the singular values sigma, which are
791 : the ratios c_i/s_i, are returned in the arguments of DSSolve().
792 : Note that the number of columns of the returned X, U, V may be smaller
793 : in the case that some c_i or s_i are zero.
794 :
795 : The number of rows of A (and U) is the value n passed with DSSetDimensions().
796 : The number of columns m and the number of rows of B (and V) must be
797 : set via DSGSVDSetDimensions().
798 :
799 : Internally, LAPACK's representation is used, U'*A*Q = C*[0 R], V'*B*Q = S*[0 R],
800 : where X = Q*inv(R) is computed at the end of DSSolve().
801 :
802 : If the compact storage format is selected, then a simplified problem is
803 : solved, where A and B are bidiagonal (possibly with an arrow), and [A;B]
804 : is assumed to have orthonormal columns. We consider two cases: (1) A and B
805 : are square mxm upper bidiagonal, and (2) A is lower bidiagonal with m+1
806 : rows and B is square upper bidiagonal. In these cases, R=I so it
807 : corresponds to the CS decomposition. The first matrix is stored in two
808 : diagonals of DS_MAT_T, while the second matrix is stored in DS_MAT_D
809 : and the remaining diagonal of DS_MAT_T.
810 :
811 : Allowed arguments of DSVectors() are DS_MAT_U, DS_MAT_V and DS_MAT_X.
812 :
813 : Used DS matrices:
814 : + DS_MAT_A - first problem matrix
815 : . DS_MAT_B - second problem matrix
816 : . DS_MAT_T - first upper bidiagonal matrix (if compact storage is selected)
817 : . DS_MAT_D - second upper bidiagonal matrix (if compact storage is selected)
818 : . DS_MAT_U - (upper) left generalized singular vectors
819 : . DS_MAT_V - (lower) left generalized singular vectors
820 : - DS_MAT_X - right generalized singular vectors
821 :
822 : Implemented methods:
823 : . 0 - Lapack (_ggsvd3 if available, or _ggsvd)
824 :
825 : .seealso: DSCreate(), DSSetType(), DSType, DSGSVDSetDimensions()
826 : M*/
827 66 : SLEPC_EXTERN PetscErrorCode DSCreate_GSVD(DS ds)
828 : {
829 66 : DS_GSVD *ctx;
830 :
831 66 : PetscFunctionBegin;
832 66 : PetscCall(PetscNew(&ctx));
833 66 : ds->data = (void*)ctx;
834 :
835 66 : ds->ops->allocate = DSAllocate_GSVD;
836 66 : ds->ops->view = DSView_GSVD;
837 66 : ds->ops->vectors = DSVectors_GSVD;
838 66 : ds->ops->sort = DSSort_GSVD;
839 66 : ds->ops->solve[0] = DSSolve_GSVD;
840 : #if !defined(PETSC_HAVE_MPIUNI)
841 66 : ds->ops->synchronize = DSSynchronize_GSVD;
842 : #endif
843 66 : ds->ops->truncate = DSTruncate_GSVD;
844 66 : ds->ops->update = DSUpdateExtraRow_GSVD;
845 66 : ds->ops->cond = DSCond_GSVD;
846 66 : ds->ops->matgetsize = DSMatGetSize_GSVD;
847 66 : ds->ops->destroy = DSDestroy_GSVD;
848 66 : ds->ops->reallocate = DSReallocate_GSVD;
849 66 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDSetDimensions_C",DSGSVDSetDimensions_GSVD));
850 66 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDGetDimensions_C",DSGSVDGetDimensions_GSVD));
851 66 : PetscFunctionReturn(PETSC_SUCCESS);
852 : }
|
