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 :
11 : #include <slepc/private/dsimpl.h> /*I "slepcds.h" I*/
12 : #include <slepcblaslapack.h>
13 :
14 : typedef struct {
15 : PetscInt m; /* number of columns */
16 : PetscInt t; /* number of rows of V after truncating */
17 : } DS_SVD;
18 :
19 89 : static PetscErrorCode DSAllocate_SVD(DS ds,PetscInt ld)
20 : {
21 89 : PetscFunctionBegin;
22 89 : if (!ds->compact) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
23 89 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
24 89 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_V));
25 89 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
26 89 : PetscCall(PetscFree(ds->perm));
27 89 : PetscCall(PetscMalloc1(ld,&ds->perm));
28 89 : PetscFunctionReturn(PETSC_SUCCESS);
29 : }
30 :
31 : /* 0 l k m-1
32 : -----------------------------------------
33 : |* . . |
34 : | * . . |
35 : | * . . |
36 : | * . . |
37 : | o o |
38 : | o o |
39 : | o o |
40 : | o o |
41 : | o o |
42 : | o o |
43 : | o x |
44 : | x x |
45 : | x x |
46 : | x x |
47 : | x x |
48 : | x x |
49 : | x x |
50 : | x x |
51 : | x x|
52 : n-1 | x|
53 : -----------------------------------------
54 : */
55 :
56 2 : static PetscErrorCode DSSwitchFormat_SVD(DS ds)
57 : {
58 2 : DS_SVD *ctx = (DS_SVD*)ds->data;
59 2 : PetscReal *T;
60 2 : PetscScalar *A;
61 2 : PetscInt i,m=ctx->m,k=ds->k,ld=ds->ld;
62 :
63 2 : PetscFunctionBegin;
64 2 : PetscCheck(m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
65 2 : PetscCheck(ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Must have compact storage");
66 : /* switch from compact (arrow) to dense storage */
67 2 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
68 2 : PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_A],&A));
69 2 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
70 2 : PetscCall(PetscArrayzero(A,ld*ld));
71 12 : for (i=0;i<k;i++) {
72 10 : A[i+i*ld] = T[i];
73 10 : A[i+k*ld] = T[i+ld];
74 : }
75 2 : A[k+k*ld] = T[k];
76 10 : for (i=k+1;i<m;i++) {
77 8 : A[i+i*ld] = T[i];
78 8 : A[i-1+i*ld] = T[i-1+ld];
79 : }
80 2 : PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_A],&A));
81 2 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
82 2 : PetscFunctionReturn(PETSC_SUCCESS);
83 : }
84 :
85 21 : static PetscErrorCode DSView_SVD(DS ds,PetscViewer viewer)
86 : {
87 21 : DS_SVD *ctx = (DS_SVD*)ds->data;
88 21 : PetscViewerFormat format;
89 21 : PetscInt i,j,r,c,m=ctx->m,rows,cols;
90 21 : PetscReal *T,value;
91 21 : const char *methodname[] = {
92 : "Implicit zero-shift QR for bidiagonals (_bdsqr)",
93 : "Divide and Conquer (_bdsdc or _gesdd)"
94 : };
95 21 : const int nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);
96 :
97 21 : PetscFunctionBegin;
98 21 : PetscCall(PetscViewerGetFormat(viewer,&format));
99 21 : if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
100 11 : PetscCall(PetscViewerASCIIPrintf(viewer,"number of columns: %" PetscInt_FMT "\n",m));
101 11 : if (ds->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer,"solving the problem with: %s\n",methodname[ds->method]));
102 11 : PetscFunctionReturn(PETSC_SUCCESS);
103 : }
104 10 : PetscCheck(m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
105 10 : if (ds->compact) {
106 6 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
107 6 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
108 6 : rows = ds->n;
109 6 : cols = ds->extrarow? m+1: m;
110 6 : if (format == PETSC_VIEWER_ASCII_MATLAB) {
111 6 : PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rows,cols));
112 6 : PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",2*ds->n));
113 6 : PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
114 66 : for (i=0;i<PetscMin(ds->n,m);i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)T[i]));
115 62 : for (i=0;i<cols-1;i++) {
116 56 : r = PetscMax(i+2,ds->k+1);
117 56 : c = i+1;
118 56 : PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",c,r,(double)T[i+ds->ld]));
119 : }
120 6 : PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_T]));
121 : } else {
122 0 : for (i=0;i<rows;i++) {
123 0 : for (j=0;j<cols;j++) {
124 0 : if (i==j) value = T[i];
125 0 : else if (i<ds->k && j==ds->k) value = T[PetscMin(i,j)+ds->ld];
126 0 : else if (i+1==j && i>=ds->k) value = T[i+ds->ld];
127 : else value = 0.0;
128 0 : PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
129 : }
130 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
131 : }
132 : }
133 6 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
134 6 : PetscCall(PetscViewerFlush(viewer));
135 6 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
136 4 : } else PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
137 10 : if (ds->state>DS_STATE_INTERMEDIATE) {
138 0 : PetscCall(DSViewMat(ds,viewer,DS_MAT_U));
139 0 : PetscCall(DSViewMat(ds,viewer,DS_MAT_V));
140 : }
141 10 : PetscFunctionReturn(PETSC_SUCCESS);
142 : }
143 :
144 4 : static PetscErrorCode DSVectors_SVD(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
145 : {
146 4 : PetscFunctionBegin;
147 4 : switch (mat) {
148 4 : case DS_MAT_U:
149 : case DS_MAT_V:
150 4 : if (rnorm) *rnorm = 0.0;
151 4 : break;
152 0 : default:
153 0 : SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
154 : }
155 4 : PetscFunctionReturn(PETSC_SUCCESS);
156 : }
157 :
158 1659 : static PetscErrorCode DSSort_SVD(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
159 : {
160 1659 : DS_SVD *ctx = (DS_SVD*)ds->data;
161 1659 : PetscInt n,l,i,*perm,ld=ds->ld;
162 1659 : PetscScalar *A;
163 1659 : PetscReal *d;
164 :
165 1659 : PetscFunctionBegin;
166 1659 : if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
167 1659 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
168 1659 : l = ds->l;
169 1659 : n = PetscMin(ds->n,ctx->m);
170 1659 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
171 1659 : perm = ds->perm;
172 1659 : if (!rr) PetscCall(DSSortEigenvaluesReal_Private(ds,d,perm));
173 0 : else PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
174 18162 : for (i=l;i<n;i++) wr[i] = d[perm[i]];
175 1659 : PetscCall(DSPermuteBoth_Private(ds,l,n,ds->n,ctx->m,DS_MAT_U,DS_MAT_V,perm));
176 18162 : for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
177 1659 : if (!ds->compact) {
178 364 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
179 4383 : for (i=l;i<n;i++) A[i+i*ld] = wr[i];
180 364 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
181 : }
182 1659 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
183 1659 : PetscFunctionReturn(PETSC_SUCCESS);
184 : }
185 :
186 1293 : static PetscErrorCode DSUpdateExtraRow_SVD(DS ds)
187 : {
188 1293 : DS_SVD *ctx = (DS_SVD*)ds->data;
189 1293 : PetscInt i;
190 1293 : PetscBLASInt n=0,m=0,ld,incx=1;
191 1293 : PetscScalar *A,*x,*y,one=1.0,zero=0.0;
192 1293 : PetscReal *T,*e,beta;
193 1293 : const PetscScalar *U;
194 :
195 1293 : PetscFunctionBegin;
196 1293 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
197 1293 : PetscCall(PetscBLASIntCast(ds->n,&n));
198 1293 : PetscCall(PetscBLASIntCast(ctx->m,&m));
199 1293 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
200 1293 : PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_U],&U));
201 1293 : if (ds->compact) {
202 1291 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
203 1291 : e = T+ld;
204 1291 : beta = e[m-1]; /* in compact, we assume all entries are zero except the last one */
205 14309 : for (i=0;i<n;i++) e[i] = PetscRealPart(beta*U[n-1+i*ld]);
206 1291 : ds->k = m;
207 1291 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
208 : } else {
209 2 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
210 2 : PetscCall(DSAllocateWork_Private(ds,2*ld,0,0));
211 2 : x = ds->work;
212 2 : y = ds->work+ld;
213 32 : for (i=0;i<n;i++) x[i] = PetscConj(A[i+m*ld]);
214 2 : PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,U,&ld,x,&incx,&zero,y,&incx));
215 32 : for (i=0;i<n;i++) A[i+m*ld] = PetscConj(y[i]);
216 2 : ds->k = m;
217 2 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
218 : }
219 1293 : PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_U],&U));
220 1293 : PetscFunctionReturn(PETSC_SUCCESS);
221 : }
222 :
223 705 : static PetscErrorCode DSTruncate_SVD(DS ds,PetscInt n,PetscBool trim)
224 : {
225 705 : PetscInt i,ld=ds->ld,l=ds->l;
226 705 : PetscScalar *A;
227 705 : DS_SVD *ctx = (DS_SVD*)ds->data;
228 :
229 705 : PetscFunctionBegin;
230 705 : if (!ds->compact && ds->extrarow) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
231 705 : if (trim) {
232 45 : if (!ds->compact && ds->extrarow) { /* clean extra column */
233 0 : for (i=l;i<ds->n;i++) A[i+ctx->m*ld] = 0.0;
234 : }
235 45 : ds->l = 0;
236 45 : ds->k = 0;
237 45 : ds->n = n;
238 45 : ctx->m = n;
239 45 : ds->t = ds->n; /* truncated length equal to the new dimension */
240 45 : ctx->t = ctx->m; /* must also keep the previous dimension of V */
241 : } else {
242 660 : if (!ds->compact && ds->extrarow && ds->k==ds->n) {
243 : /* copy entries of extra column to the new position, then clean last row */
244 0 : for (i=l;i<n;i++) A[i+n*ld] = A[i+ctx->m*ld];
245 0 : for (i=l;i<ds->n;i++) A[i+ctx->m*ld] = 0.0;
246 : }
247 660 : ds->k = (ds->extrarow)? n: 0;
248 660 : ds->t = ds->n; /* truncated length equal to previous dimension */
249 660 : ctx->t = ctx->m; /* must also keep the previous dimension of V */
250 660 : ds->n = n;
251 660 : ctx->m = n;
252 : }
253 705 : if (!ds->compact && ds->extrarow) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
254 705 : PetscFunctionReturn(PETSC_SUCCESS);
255 : }
256 :
257 : /*
258 : DSArrowBidiag reduces a real square arrowhead matrix of the form
259 :
260 : [ d 0 0 0 e ]
261 : [ 0 d 0 0 e ]
262 : A = [ 0 0 d 0 e ]
263 : [ 0 0 0 d e ]
264 : [ 0 0 0 0 d ]
265 :
266 : to upper bidiagonal form
267 :
268 : [ d e 0 0 0 ]
269 : [ 0 d e 0 0 ]
270 : B = Q'*A*P = [ 0 0 d e 0 ],
271 : [ 0 0 0 d e ]
272 : [ 0 0 0 0 d ]
273 :
274 : where P,Q are orthogonal matrices. Uses plane rotations with a bulge chasing scheme.
275 : On input, P and Q must be initialized to the identity matrix.
276 : */
277 662 : static PetscErrorCode DSArrowBidiag(PetscBLASInt n,PetscReal *d,PetscReal *e,PetscScalar *Q,PetscBLASInt ldq,PetscScalar *P,PetscBLASInt ldp)
278 : {
279 662 : PetscBLASInt i,j,j2,one=1;
280 662 : PetscReal c,s,ct,st,off,temp0,temp1,temp2;
281 :
282 662 : PetscFunctionBegin;
283 662 : if (n<=2) PetscFunctionReturn(PETSC_SUCCESS);
284 :
285 3095 : for (j=0;j<n-2;j++) {
286 :
287 : /* Eliminate entry e(j) by a rotation in the planes (j,j+1) */
288 2433 : temp0 = e[j+1];
289 2433 : PetscCallBLAS("LAPACKlartg",LAPACKREALlartg_(&temp0,&e[j],&c,&s,&e[j+1]));
290 2433 : s = -s;
291 :
292 : /* Apply rotation to Q */
293 2433 : j2 = j+2;
294 2433 : PetscCallBLAS("BLASrot",BLASMIXEDrot_(&j2,Q+j*ldq,&one,Q+(j+1)*ldq,&one,&c,&s));
295 :
296 : /* Apply rotation to diagonal elements, eliminate newly introduced entry A(j+1,j) */
297 2433 : temp0 = d[j+1];
298 2433 : temp1 = c*temp0;
299 2433 : temp2 = -s*d[j];
300 2433 : PetscCallBLAS("LAPACKlartg",LAPACKREALlartg_(&temp1,&temp2,&ct,&st,&d[j+1]));
301 2433 : st = -st;
302 2433 : e[j] = -c*st*d[j] + s*ct*temp0;
303 2433 : d[j] = c*ct*d[j] + s*st*temp0;
304 :
305 : /* Apply rotation to P */
306 2433 : PetscCallBLAS("BLASrot",BLASMIXEDrot_(&j2,P+j*ldp,&one,P+(j+1)*ldp,&one,&ct,&st));
307 :
308 : /* Chase newly introduced off-diagonal entry to the top left corner */
309 5928 : for (i=j-1;i>=0;i--) {
310 :
311 : /* Upper bulge */
312 3495 : off = -st*e[i];
313 3495 : e[i] = ct*e[i];
314 3495 : temp0 = e[i+1];
315 3495 : PetscCallBLAS("LAPACKlartg",LAPACKREALlartg_(&temp0,&off,&c,&s,&e[i+1]));
316 3495 : s = -s;
317 3495 : PetscCallBLAS("BLASrot",BLASMIXEDrot_(&j2,Q+i*ldq,&one,Q+(i+1)*ldq,&one,&c,&s));
318 :
319 : /* Lower bulge */
320 3495 : temp0 = d[i+1];
321 3495 : temp1 = -s*e[i] + c*temp0;
322 3495 : temp2 = c*e[i] + s*temp0;
323 3495 : off = -s*d[i];
324 3495 : PetscCallBLAS("LAPACKlartg",LAPACKREALlartg_(&temp1,&off,&ct,&st,&d[i+1]));
325 3495 : st = -st;
326 3495 : e[i] = -c*st*d[i] + ct*temp2;
327 3495 : d[i] = c*ct*d[i] + st*temp2;
328 3495 : PetscCallBLAS("BLASrot",BLASMIXEDrot_(&j2,P+i*ldp,&one,P+(i+1)*ldp,&one,&ct,&st));
329 : }
330 : }
331 662 : PetscFunctionReturn(PETSC_SUCCESS);
332 : }
333 :
334 : /*
335 : Reduce to bidiagonal form by means of DSArrowBidiag.
336 : */
337 1706 : static PetscErrorCode DSIntermediate_SVD(DS ds)
338 : {
339 1706 : DS_SVD *ctx = (DS_SVD*)ds->data;
340 1706 : PetscInt i,j;
341 1706 : PetscBLASInt n1 = 0,n2,m2,lwork,info,l = 0,n = 0,m = 0,nm,ld,off;
342 1706 : PetscScalar *A,*U,*V,*W,*work,*tauq,*taup;
343 1706 : PetscReal *d,*e;
344 :
345 1706 : PetscFunctionBegin;
346 1706 : PetscCall(PetscBLASIntCast(ds->n,&n));
347 1706 : PetscCall(PetscBLASIntCast(ctx->m,&m));
348 1706 : PetscCall(PetscBLASIntCast(ds->l,&l));
349 1706 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
350 1706 : PetscCall(PetscBLASIntCast(PetscMax(0,ds->k-l+1),&n1)); /* size of leading block, excl. locked */
351 1706 : n2 = n-l; /* n2 = n1 + size of trailing block */
352 1706 : m2 = m-l;
353 1706 : off = l+l*ld;
354 1706 : nm = PetscMin(n,m);
355 1706 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
356 1706 : e = d+ld;
357 1706 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
358 1706 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
359 1706 : PetscCall(PetscArrayzero(U,ld*ld));
360 19801 : for (i=0;i<n;i++) U[i+i*ld] = 1.0;
361 1706 : PetscCall(PetscArrayzero(V,ld*ld));
362 19773 : for (i=0;i<m;i++) V[i+i*ld] = 1.0;
363 :
364 1706 : if (ds->compact) {
365 :
366 1292 : if (ds->state<DS_STATE_INTERMEDIATE) PetscCall(DSArrowBidiag(n1,d+l,e+l,U+off,ld,V+off,ld));
367 :
368 : } else {
369 :
370 414 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
371 565 : for (i=0;i<l;i++) { d[i] = PetscRealPart(A[i+i*ld]); e[i] = 0.0; }
372 :
373 414 : if (ds->state<DS_STATE_INTERMEDIATE) {
374 414 : lwork = (m+n)*16;
375 414 : PetscCall(DSAllocateWork_Private(ds,2*nm+ld*ld+lwork,0,0));
376 414 : tauq = ds->work;
377 414 : taup = ds->work+nm;
378 414 : W = ds->work+2*nm;
379 414 : work = ds->work+2*nm+ld*ld;
380 5453 : for (j=0;j<m;j++) PetscCall(PetscArraycpy(W+j*ld,A+j*ld,n));
381 414 : PetscCallBLAS("LAPACKgebrd",LAPACKgebrd_(&n2,&m2,W+off,&ld,d+l,e+l,tauq,taup,work,&lwork,&info));
382 414 : SlepcCheckLapackInfo("gebrd",info);
383 414 : PetscCallBLAS("LAPACKormbr",LAPACKormbr_("Q","L","N",&n2,&n2,&m2,W+off,&ld,tauq,U+off,&ld,work,&lwork,&info));
384 414 : SlepcCheckLapackInfo("ormbr",info);
385 414 : PetscCallBLAS("LAPACKormbr",LAPACKormbr_("P","R","N",&m2,&m2,&n2,W+off,&ld,taup,V+off,&ld,work,&lwork,&info));
386 414 : SlepcCheckLapackInfo("ormbr",info);
387 : } else {
388 : /* copy bidiagonal to d,e */
389 0 : for (i=l;i<nm;i++) d[i] = PetscRealPart(A[i+i*ld]);
390 0 : for (i=l;i<nm-1;i++) e[i] = PetscRealPart(A[i+(i+1)*ld]);
391 : }
392 414 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
393 : }
394 1706 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
395 1706 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
396 1706 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
397 1706 : PetscFunctionReturn(PETSC_SUCCESS);
398 : }
399 :
400 1706 : static PetscErrorCode DSSolve_SVD_QR(DS ds,PetscScalar *wr,PetscScalar *wi)
401 : {
402 1706 : DS_SVD *ctx = (DS_SVD*)ds->data;
403 1706 : PetscInt i,j;
404 1706 : PetscBLASInt n1,m1,info,l = 0,n = 0,m = 0,nm,ld,off,zero=0;
405 1706 : PetscScalar *A,*U,*V,*Vt;
406 1706 : PetscReal *d,*e;
407 :
408 1706 : PetscFunctionBegin;
409 1706 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
410 1706 : PetscCall(PetscBLASIntCast(ds->n,&n));
411 1706 : PetscCall(PetscBLASIntCast(ctx->m,&m));
412 1706 : PetscCall(PetscBLASIntCast(ds->l,&l));
413 1706 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
414 1706 : n1 = n-l; /* n1 = size of leading block, excl. locked + size of trailing block */
415 1706 : m1 = m-l;
416 1706 : nm = PetscMin(n1,m1);
417 1706 : off = l+l*ld;
418 1706 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
419 1706 : e = d+ld;
420 :
421 : /* Reduce to bidiagonal form */
422 1706 : PetscCall(DSIntermediate_SVD(ds));
423 :
424 1706 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
425 1706 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
426 :
427 : /* solve bidiagonal SVD problem */
428 2427 : for (i=0;i<l;i++) wr[i] = d[i];
429 1706 : PetscCall(DSAllocateWork_Private(ds,ld*ld,4*n1,0));
430 1706 : Vt = ds->work;
431 19052 : for (i=l;i<m;i++) {
432 209724 : for (j=l;j<m;j++) {
433 192378 : Vt[i+j*ld] = PetscConj(V[j+i*ld]); /* Lapack expects transposed VT */
434 : }
435 : }
436 1709 : PetscCallBLAS("LAPACKbdsqr",LAPACKbdsqr_(n>=m?"U":"L",&nm,&m1,&n1,&zero,d+l,e+l,Vt+off,&ld,U+off,&ld,NULL,&ld,ds->rwork,&info));
437 1706 : SlepcCheckLapackInfo("bdsqr",info);
438 19052 : for (i=l;i<m;i++) {
439 209724 : for (j=l;j<m;j++) {
440 192378 : V[i+j*ld] = PetscConj(Vt[j+i*ld]); /* transpose VT returned by Lapack */
441 : }
442 : }
443 19046 : for (i=l;i<PetscMin(ds->n,ctx->m);i++) wr[i] = d[i];
444 :
445 : /* create diagonal matrix as a result */
446 1706 : if (ds->compact) PetscCall(PetscArrayzero(e,n-1));
447 : else {
448 414 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
449 5302 : for (i=l;i<m;i++) PetscCall(PetscArrayzero(A+l+i*ld,n-l));
450 5296 : for (i=l;i<PetscMin(ds->n,ctx->m);i++) A[i+i*ld] = d[i];
451 414 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
452 : }
453 1706 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
454 1706 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
455 1706 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
456 1706 : PetscFunctionReturn(PETSC_SUCCESS);
457 : }
458 :
459 5 : static PetscErrorCode DSSolve_SVD_DC(DS ds,PetscScalar *wr,PetscScalar *wi)
460 : {
461 5 : DS_SVD *ctx = (DS_SVD*)ds->data;
462 5 : PetscInt i,j;
463 5 : PetscBLASInt n1,m1,info,l = 0,n = 0,m = 0,nm,ld,off,lwork;
464 5 : PetscScalar *A,*U,*V,*W,qwork;
465 5 : PetscReal *d,*e,*Ur,*Vr;
466 :
467 5 : PetscFunctionBegin;
468 5 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
469 5 : PetscCall(PetscBLASIntCast(ds->n,&n));
470 5 : PetscCall(PetscBLASIntCast(ctx->m,&m));
471 5 : PetscCall(PetscBLASIntCast(ds->l,&l));
472 5 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
473 5 : n1 = n-l; /* n1 = size of leading block, excl. locked + size of trailing block */
474 5 : m1 = m-l;
475 5 : off = l+l*ld;
476 5 : if (ds->compact) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
477 5 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
478 5 : PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_U],&U));
479 5 : PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_V],&V));
480 5 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
481 5 : e = d+ld;
482 5 : PetscCall(PetscArrayzero(U,ld*ld));
483 9 : for (i=0;i<l;i++) U[i+i*ld] = 1.0;
484 5 : PetscCall(PetscArrayzero(V,ld*ld));
485 9 : for (i=0;i<l;i++) V[i+i*ld] = 1.0;
486 :
487 5 : if (ds->state>DS_STATE_RAW) {
488 : /* solve bidiagonal SVD problem */
489 1 : for (i=0;i<l;i++) wr[i] = d[i];
490 : #if defined(PETSC_USE_COMPLEX)
491 : PetscCall(DSAllocateWork_Private(ds,0,3*n1*n1+4*n1+2*ld*ld,8*n1));
492 : Ur = ds->rwork+3*n1*n1+4*n1;
493 : Vr = ds->rwork+3*n1*n1+4*n1+ld*ld;
494 : #else
495 1 : PetscCall(DSAllocateWork_Private(ds,0,3*n1*n1+4*n1+ld*ld,8*n1));
496 1 : Ur = U;
497 1 : Vr = ds->rwork+3*n1*n1+4*n1;
498 : #endif
499 1 : PetscCallBLAS("LAPACKbdsdc",LAPACKbdsdc_("U","I",&n1,d+l,e+l,Ur+off,&ld,Vr+off,&ld,NULL,NULL,ds->rwork,ds->iwork,&info));
500 1 : SlepcCheckLapackInfo("bdsdc",info);
501 11 : for (i=l;i<n;i++) {
502 110 : for (j=l;j<n;j++) {
503 : #if defined(PETSC_USE_COMPLEX)
504 : U[i+j*ld] = Ur[i+j*ld];
505 : #endif
506 100 : V[i+j*ld] = PetscConj(Vr[j+i*ld]); /* transpose VT returned by Lapack */
507 : }
508 : }
509 : } else {
510 : /* solve general rectangular SVD problem */
511 4 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
512 4 : PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&W));
513 4 : if (ds->compact) PetscCall(DSSwitchFormat_SVD(ds));
514 8 : for (i=0;i<l;i++) wr[i] = d[i];
515 4 : nm = PetscMin(n,m);
516 4 : PetscCall(DSAllocateWork_Private(ds,0,0,8*nm));
517 4 : lwork = -1;
518 : #if defined(PETSC_USE_COMPLEX)
519 : PetscCall(DSAllocateWork_Private(ds,0,5*nm*nm+7*nm,0));
520 : PetscCallBLAS("LAPACKgesdd",LAPACKgesdd_("A",&n1,&m1,A+off,&ld,d+l,U+off,&ld,W+off,&ld,&qwork,&lwork,ds->rwork,ds->iwork,&info));
521 : #else
522 4 : PetscCallBLAS("LAPACKgesdd",LAPACKgesdd_("A",&n1,&m1,A+off,&ld,d+l,U+off,&ld,W+off,&ld,&qwork,&lwork,ds->iwork,&info));
523 : #endif
524 4 : SlepcCheckLapackInfo("gesdd",info);
525 4 : PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(qwork),&lwork));
526 4 : PetscCall(DSAllocateWork_Private(ds,lwork,0,0));
527 : #if defined(PETSC_USE_COMPLEX)
528 : PetscCallBLAS("LAPACKgesdd",LAPACKgesdd_("A",&n1,&m1,A+off,&ld,d+l,U+off,&ld,W+off,&ld,ds->work,&lwork,ds->rwork,ds->iwork,&info));
529 : #else
530 4 : PetscCallBLAS("LAPACKgesdd",LAPACKgesdd_("A",&n1,&m1,A+off,&ld,d+l,U+off,&ld,W+off,&ld,ds->work,&lwork,ds->iwork,&info));
531 : #endif
532 4 : SlepcCheckLapackInfo("gesdd",info);
533 40 : for (i=l;i<m;i++) {
534 364 : for (j=l;j<m;j++) V[i+j*ld] = PetscConj(W[j+i*ld]); /* transpose VT returned by Lapack */
535 : }
536 4 : PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&W));
537 : }
538 51 : for (i=l;i<PetscMin(ds->n,ctx->m);i++) wr[i] = d[i];
539 :
540 : /* create diagonal matrix as a result */
541 5 : if (ds->compact) PetscCall(PetscArrayzero(e,n-1));
542 : else {
543 22 : for (i=l;i<m;i++) PetscCall(PetscArrayzero(A+l+i*ld,n-l));
544 32 : for (i=l;i<n;i++) A[i+i*ld] = d[i];
545 : }
546 5 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
547 5 : PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_U],&U));
548 5 : PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_V],&V));
549 5 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
550 5 : PetscFunctionReturn(PETSC_SUCCESS);
551 : }
552 :
553 : #if !defined(PETSC_HAVE_MPIUNI)
554 6 : static PetscErrorCode DSSynchronize_SVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
555 : {
556 6 : PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
557 6 : PetscMPIInt n,rank,off=0,size,ldn,ld3;
558 6 : PetscScalar *A,*U,*V;
559 6 : PetscReal *T;
560 :
561 6 : PetscFunctionBegin;
562 6 : if (ds->compact) kr = 3*ld;
563 0 : else k = (ds->n-l)*ld;
564 6 : if (ds->state>DS_STATE_RAW) k += 2*(ds->n-l)*ld;
565 6 : if (eigr) k += ds->n-l;
566 6 : PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
567 6 : PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
568 6 : PetscCall(PetscMPIIntCast(ds->n-l,&n));
569 6 : PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
570 6 : PetscCall(PetscMPIIntCast(3*ld,&ld3));
571 6 : if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
572 0 : else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
573 6 : if (ds->state>DS_STATE_RAW) {
574 6 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
575 6 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
576 : }
577 6 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
578 6 : if (!rank) {
579 3 : if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
580 0 : else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
581 3 : if (ds->state>DS_STATE_RAW) {
582 3 : PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
583 3 : PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
584 : }
585 3 : if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
586 : }
587 12 : PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
588 6 : if (rank) {
589 3 : if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
590 0 : else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
591 3 : if (ds->state>DS_STATE_RAW) {
592 3 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
593 3 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
594 : }
595 3 : if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
596 : }
597 6 : if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
598 0 : else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
599 6 : if (ds->state>DS_STATE_RAW) {
600 6 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
601 6 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
602 : }
603 6 : PetscFunctionReturn(PETSC_SUCCESS);
604 : }
605 : #endif
606 :
607 3680 : static PetscErrorCode DSMatGetSize_SVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
608 : {
609 3680 : DS_SVD *ctx = (DS_SVD*)ds->data;
610 :
611 3680 : PetscFunctionBegin;
612 3680 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
613 3680 : switch (t) {
614 416 : case DS_MAT_A:
615 416 : *rows = ds->n;
616 416 : *cols = ds->extrarow? ctx->m+1: ctx->m;
617 416 : break;
618 0 : case DS_MAT_T:
619 0 : *rows = ds->n;
620 0 : *cols = PetscDefined(USE_COMPLEX)? 2: 3;
621 0 : break;
622 1631 : case DS_MAT_U:
623 1631 : *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
624 1631 : *cols = ds->n;
625 1631 : break;
626 1633 : case DS_MAT_V:
627 1633 : *rows = ds->state==DS_STATE_TRUNCATED? ctx->t: ctx->m;
628 1633 : *cols = ctx->m;
629 1633 : break;
630 0 : default:
631 0 : SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
632 : }
633 3680 : PetscFunctionReturn(PETSC_SUCCESS);
634 : }
635 :
636 1711 : static PetscErrorCode DSSVDSetDimensions_SVD(DS ds,PetscInt m)
637 : {
638 1711 : DS_SVD *ctx = (DS_SVD*)ds->data;
639 :
640 1711 : PetscFunctionBegin;
641 1711 : DSCheckAlloc(ds,1);
642 1711 : if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
643 0 : ctx->m = ds->ld;
644 : } else {
645 1711 : PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
646 1711 : ctx->m = m;
647 : }
648 1711 : PetscFunctionReturn(PETSC_SUCCESS);
649 : }
650 :
651 : /*@
652 : DSSVDSetDimensions - Sets the number of columns for a DSSVD.
653 :
654 : Logically Collective
655 :
656 : Input Parameters:
657 : + ds - the direct solver context
658 : - m - the number of columns
659 :
660 : Notes:
661 : This call is complementary to DSSetDimensions(), to provide a dimension
662 : that is specific to this DS type.
663 :
664 : Level: intermediate
665 :
666 : .seealso: DSSVDGetDimensions(), DSSetDimensions()
667 : @*/
668 1711 : PetscErrorCode DSSVDSetDimensions(DS ds,PetscInt m)
669 : {
670 1711 : PetscFunctionBegin;
671 1711 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
672 5133 : PetscValidLogicalCollectiveInt(ds,m,2);
673 1711 : PetscTryMethod(ds,"DSSVDSetDimensions_C",(DS,PetscInt),(ds,m));
674 1711 : PetscFunctionReturn(PETSC_SUCCESS);
675 : }
676 :
677 4 : static PetscErrorCode DSSVDGetDimensions_SVD(DS ds,PetscInt *m)
678 : {
679 4 : DS_SVD *ctx = (DS_SVD*)ds->data;
680 :
681 4 : PetscFunctionBegin;
682 4 : *m = ctx->m;
683 4 : PetscFunctionReturn(PETSC_SUCCESS);
684 : }
685 :
686 : /*@
687 : DSSVDGetDimensions - Returns the number of columns for a DSSVD.
688 :
689 : Not Collective
690 :
691 : Input Parameter:
692 : . ds - the direct solver context
693 :
694 : Output Parameters:
695 : . m - the number of columns
696 :
697 : Level: intermediate
698 :
699 : .seealso: DSSVDSetDimensions()
700 : @*/
701 4 : PetscErrorCode DSSVDGetDimensions(DS ds,PetscInt *m)
702 : {
703 4 : PetscFunctionBegin;
704 4 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
705 4 : PetscAssertPointer(m,2);
706 4 : PetscUseMethod(ds,"DSSVDGetDimensions_C",(DS,PetscInt*),(ds,m));
707 4 : PetscFunctionReturn(PETSC_SUCCESS);
708 : }
709 :
710 86 : static PetscErrorCode DSDestroy_SVD(DS ds)
711 : {
712 86 : PetscFunctionBegin;
713 86 : PetscCall(PetscFree(ds->data));
714 86 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSSVDSetDimensions_C",NULL));
715 86 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSSVDGetDimensions_C",NULL));
716 86 : PetscFunctionReturn(PETSC_SUCCESS);
717 : }
718 :
719 6 : static PetscErrorCode DSSetCompact_SVD(DS ds,PetscBool comp)
720 : {
721 6 : PetscFunctionBegin;
722 6 : if (!comp) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
723 6 : PetscFunctionReturn(PETSC_SUCCESS);
724 : }
725 :
726 : /*MC
727 : DSSVD - Dense Singular Value Decomposition.
728 :
729 : Level: beginner
730 :
731 : Notes:
732 : The problem is expressed as A = U*Sigma*V', where A is rectangular in
733 : general, with n rows and m columns. Sigma is a diagonal matrix whose diagonal
734 : elements are the arguments of DSSolve(). After solve, A is overwritten
735 : with Sigma.
736 :
737 : The orthogonal (or unitary) matrices of left and right singular vectors, U
738 : and V, have size n and m, respectively. The number of columns m must
739 : be specified via DSSVDSetDimensions().
740 :
741 : If the DS object is in the intermediate state, A is assumed to be in upper
742 : bidiagonal form (possibly with an arrow) and is stored in compact format
743 : on matrix T. Otherwise, no particular structure is assumed. The compact
744 : storage is implemented for the square case only, m=n. The extra row should
745 : be interpreted in this case as an extra column.
746 :
747 : Used DS matrices:
748 : + DS_MAT_A - problem matrix (used only if compact=false)
749 : . DS_MAT_T - upper bidiagonal matrix
750 : . DS_MAT_U - left singular vectors
751 : - DS_MAT_V - right singular vectors
752 :
753 : Implemented methods:
754 : + 0 - Implicit zero-shift QR for bidiagonals (_bdsqr)
755 : - 1 - Divide and Conquer (_bdsdc or _gesdd)
756 :
757 : .seealso: DSCreate(), DSSetType(), DSType, DSSVDSetDimensions()
758 : M*/
759 86 : SLEPC_EXTERN PetscErrorCode DSCreate_SVD(DS ds)
760 : {
761 86 : DS_SVD *ctx;
762 :
763 86 : PetscFunctionBegin;
764 86 : PetscCall(PetscNew(&ctx));
765 86 : ds->data = (void*)ctx;
766 :
767 86 : ds->ops->allocate = DSAllocate_SVD;
768 86 : ds->ops->view = DSView_SVD;
769 86 : ds->ops->vectors = DSVectors_SVD;
770 86 : ds->ops->solve[0] = DSSolve_SVD_QR;
771 86 : ds->ops->solve[1] = DSSolve_SVD_DC;
772 86 : ds->ops->sort = DSSort_SVD;
773 86 : ds->ops->truncate = DSTruncate_SVD;
774 86 : ds->ops->update = DSUpdateExtraRow_SVD;
775 86 : ds->ops->destroy = DSDestroy_SVD;
776 86 : ds->ops->matgetsize = DSMatGetSize_SVD;
777 : #if !defined(PETSC_HAVE_MPIUNI)
778 86 : ds->ops->synchronize = DSSynchronize_SVD;
779 : #endif
780 86 : ds->ops->setcompact = DSSetCompact_SVD;
781 86 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSSVDSetDimensions_C",DSSVDSetDimensions_SVD));
782 86 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSSVDGetDimensions_C",DSSVDGetDimensions_SVD));
783 86 : PetscFunctionReturn(PETSC_SUCCESS);
784 : }
|
