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 : PetscBool reorth; /* reorthogonalize left vectors */
18 : } DS_HSVD;
19 :
20 21 : static PetscErrorCode DSAllocate_HSVD(DS ds,PetscInt ld)
21 : {
22 21 : PetscFunctionBegin;
23 21 : if (!ds->compact) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
24 21 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
25 21 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_V));
26 21 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
27 21 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_D));
28 21 : PetscCall(PetscFree(ds->perm));
29 21 : PetscCall(PetscMalloc1(ld,&ds->perm));
30 21 : PetscFunctionReturn(PETSC_SUCCESS);
31 : }
32 :
33 : /* 0 l k m-1
34 : -----------------------------------------
35 : |* . . |
36 : | * . . |
37 : | * . . |
38 : | * . . |
39 : | o o |
40 : | o o |
41 : | o o |
42 : | o o |
43 : | o o |
44 : | o o |
45 : | o x |
46 : | x x |
47 : | x x |
48 : | x x |
49 : | x x |
50 : | x x |
51 : | x x |
52 : | x x |
53 : | x x|
54 : n-1 | x|
55 : -----------------------------------------
56 : */
57 :
58 10 : static PetscErrorCode DSView_HSVD(DS ds,PetscViewer viewer)
59 : {
60 10 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
61 10 : PetscViewerFormat format;
62 10 : PetscInt i,j,r,c,m=ctx->m,rows,cols;
63 10 : PetscReal *T,*S,value;
64 10 : const char *methodname[] = {
65 : "Cross product A'*Omega*A"
66 : };
67 10 : const int nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);
68 :
69 10 : PetscFunctionBegin;
70 10 : PetscCall(PetscViewerGetFormat(viewer,&format));
71 10 : if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
72 5 : if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscCall(PetscViewerASCIIPrintf(viewer,"number of columns: %" PetscInt_FMT "\n",m));
73 5 : if (ds->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer,"solving the problem with: %s\n",methodname[ds->method]));
74 5 : if (ctx->reorth) PetscCall(PetscViewerASCIIPrintf(viewer,"reorthogonalizing left vectors\n"));
75 5 : PetscFunctionReturn(PETSC_SUCCESS);
76 : }
77 5 : PetscCheck(m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
78 5 : if (ds->compact) {
79 4 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
80 4 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
81 4 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
82 4 : rows = ds->n;
83 4 : cols = ds->extrarow? m+1: m;
84 4 : if (format == PETSC_VIEWER_ASCII_MATLAB) {
85 4 : PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rows,cols));
86 4 : PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",2*ds->n));
87 4 : PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
88 44 : 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]));
89 42 : for (i=0;i<cols-1;i++) {
90 38 : c = PetscMax(i+2,ds->k+1);
91 38 : r = i+1;
92 38 : value = i<ds->l? 0.0: T[i+ds->ld];
93 38 : PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",r,c,(double)value));
94 : }
95 4 : PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_T]));
96 4 : PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",ds->n,ds->n));
97 4 : PetscCall(PetscViewerASCIIPrintf(viewer,"omega = zeros(%" PetscInt_FMT ",3);\n",3*ds->n));
98 4 : PetscCall(PetscViewerASCIIPrintf(viewer,"omega = [\n"));
99 44 : for (i=0;i<ds->n;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)S[i]));
100 4 : PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(omega);\n",DSMatName[DS_MAT_D]));
101 : } else {
102 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"T\n"));
103 0 : for (i=0;i<rows;i++) {
104 0 : for (j=0;j<cols;j++) {
105 0 : if (i==j) value = T[i];
106 0 : else if (i<ds->l) value = 0.0;
107 0 : else if (i<ds->k && j==ds->k) value = T[PetscMin(i,j)+ds->ld];
108 0 : else if (i+1==j && i>=ds->k) value = T[i+ds->ld];
109 : else value = 0.0;
110 0 : PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
111 : }
112 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
113 : }
114 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"omega\n"));
115 0 : for (i=0;i<ds->n;i++) {
116 0 : for (j=0;j<ds->n;j++) {
117 0 : if (i==j) value = S[i];
118 : else value = 0.0;
119 0 : PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
120 : }
121 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
122 : }
123 : }
124 4 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
125 4 : PetscCall(PetscViewerFlush(viewer));
126 4 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
127 4 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
128 : } else {
129 1 : PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
130 1 : PetscCall(DSViewMat(ds,viewer,DS_MAT_D));
131 : }
132 5 : if (ds->state>DS_STATE_INTERMEDIATE) {
133 0 : PetscCall(DSViewMat(ds,viewer,DS_MAT_U));
134 0 : PetscCall(DSViewMat(ds,viewer,DS_MAT_V));
135 : }
136 5 : PetscFunctionReturn(PETSC_SUCCESS);
137 : }
138 :
139 1 : static PetscErrorCode DSVectors_HSVD(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
140 : {
141 1 : PetscFunctionBegin;
142 1 : switch (mat) {
143 1 : case DS_MAT_U:
144 : case DS_MAT_V:
145 1 : if (rnorm) *rnorm = 0.0;
146 1 : break;
147 0 : default:
148 0 : SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
149 : }
150 1 : PetscFunctionReturn(PETSC_SUCCESS);
151 : }
152 :
153 383 : static PetscErrorCode DSSort_HSVD(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
154 : {
155 383 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
156 383 : PetscInt n,l,i,*perm,ld=ds->ld;
157 383 : PetscScalar *A;
158 383 : PetscReal *d,*s;
159 :
160 383 : PetscFunctionBegin;
161 383 : if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
162 383 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
163 383 : l = ds->l;
164 383 : n = PetscMin(ds->n,ctx->m);
165 383 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
166 383 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
167 383 : PetscCall(DSAllocateWork_Private(ds,0,ds->ld,0));
168 383 : perm = ds->perm;
169 383 : if (!rr) PetscCall(DSSortEigenvaluesReal_Private(ds,d,perm));
170 0 : else PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
171 383 : PetscCall(PetscArraycpy(ds->rwork,s,n));
172 10275 : for (i=l;i<n;i++) s[i] = ds->rwork[perm[i]];
173 10275 : for (i=l;i<n;i++) wr[i] = d[perm[i]];
174 383 : PetscCall(DSPermuteBoth_Private(ds,l,n,ds->n,ctx->m,DS_MAT_U,DS_MAT_V,perm));
175 10275 : for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
176 383 : if (!ds->compact) {
177 3 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
178 533 : for (i=l;i<n;i++) A[i+i*ld] = wr[i];
179 3 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
180 : }
181 383 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
182 383 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
183 383 : PetscFunctionReturn(PETSC_SUCCESS);
184 : }
185 :
186 378 : static PetscErrorCode DSUpdateExtraRow_HSVD(DS ds)
187 : {
188 378 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
189 378 : PetscInt i;
190 378 : PetscBLASInt n=0,m=0,ld,l;
191 378 : const PetscScalar *U;
192 378 : PetscReal *T,*e,*Omega,beta;
193 :
194 378 : PetscFunctionBegin;
195 378 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
196 378 : PetscCall(PetscBLASIntCast(ds->n,&n));
197 378 : PetscCall(PetscBLASIntCast(ctx->m,&m));
198 378 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
199 378 : PetscCall(PetscBLASIntCast(ds->l,&l));
200 378 : PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_U],&U));
201 378 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&Omega));
202 378 : PetscCheck(ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-compact storage");
203 378 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
204 378 : e = T+ld;
205 378 : beta = PetscAbs(e[m-1]); /* in compact, we assume all entries are zero except the last one */
206 10068 : for (i=0;i<n;i++) e[i] = PetscRealPart(beta*U[n-1+i*ld]*Omega[i]);
207 378 : ds->k = m;
208 378 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
209 378 : PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_U],&U));
210 378 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&Omega));
211 378 : PetscFunctionReturn(PETSC_SUCCESS);
212 : }
213 :
214 376 : static PetscErrorCode DSTruncate_HSVD(DS ds,PetscInt n,PetscBool trim)
215 : {
216 376 : PetscInt i,ld=ds->ld,l=ds->l;
217 376 : PetscScalar *A;
218 376 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
219 :
220 376 : PetscFunctionBegin;
221 376 : if (!ds->compact && ds->extrarow) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
222 376 : if (trim) {
223 14 : if (!ds->compact && ds->extrarow) { /* clean extra column */
224 0 : for (i=l;i<ds->n;i++) A[i+ctx->m*ld] = 0.0;
225 : }
226 14 : ds->l = 0;
227 14 : ds->k = 0;
228 14 : ds->n = n;
229 14 : ctx->m = n;
230 14 : ds->t = ds->n; /* truncated length equal to the new dimension */
231 14 : ctx->t = ctx->m; /* must also keep the previous dimension of V */
232 : } else {
233 362 : if (!ds->compact && ds->extrarow && ds->k==ds->n) {
234 : /* copy entries of extra column to the new position, then clean last row */
235 0 : for (i=l;i<n;i++) A[i+n*ld] = A[i+ctx->m*ld];
236 0 : for (i=l;i<ds->n;i++) A[i+ctx->m*ld] = 0.0;
237 : }
238 362 : ds->k = (ds->extrarow)? n: 0;
239 362 : ds->t = ds->n; /* truncated length equal to previous dimension */
240 362 : ctx->t = ctx->m; /* must also keep the previous dimension of V */
241 362 : ds->n = n;
242 362 : ctx->m = n;
243 : }
244 376 : if (!ds->compact && ds->extrarow) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
245 376 : PetscFunctionReturn(PETSC_SUCCESS);
246 : }
247 :
248 383 : static PetscErrorCode DSSolve_HSVD_CROSS(DS ds,PetscScalar *wr,PetscScalar *wi)
249 : {
250 383 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
251 383 : PetscInt i,j,k=ds->k,rwu=0,iwu=0,swu=0,nv;
252 383 : PetscBLASInt n1,n2,info,l=0,n=0,m=0,ld,off,one=1,*perm,*cmplx,incx=1,lwork;
253 383 : PetscScalar *A,*U,*V,scal,*R,sone=1.0,szero=0.0;
254 383 : PetscReal *d,*e,*dd,*ee,*Omega;
255 :
256 383 : PetscFunctionBegin;
257 383 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
258 383 : PetscCall(PetscBLASIntCast(ds->n,&n));
259 383 : PetscCall(PetscBLASIntCast(ctx->m,&m));
260 383 : PetscCheck(!ds->compact || n==m,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-square matrices in compact storage");
261 383 : PetscCheck(ds->compact || n>=m,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for the case of more columns than rows");
262 383 : PetscCall(PetscBLASIntCast(ds->l,&l));
263 383 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
264 383 : PetscCall(PetscBLASIntCast(PetscMax(0,ds->k-ds->l+1),&n2));
265 383 : n1 = n-l; /* n1 = size of leading block, excl. locked + size of trailing block */
266 383 : off = l+l*ld;
267 383 : if (!ds->compact) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
268 383 : PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_U],&U));
269 383 : PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_V],&V));
270 383 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
271 383 : e = d+ld;
272 383 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&Omega));
273 383 : PetscCall(PetscArrayzero(U,ld*ld));
274 731 : for (i=0;i<l;i++) U[i+i*ld] = 1.0;
275 383 : PetscCall(PetscArrayzero(V,ld*ld));
276 11341 : for (i=0;i<n;i++) V[i+i*ld] = 1.0;
277 731 : for (i=0;i<l;i++) wr[i] = d[i];
278 383 : if (wi) for (i=0;i<l;i++) wi[i] = 0.0;
279 :
280 383 : if (ds->compact) {
281 : /* Form the arrow tridiagonal cross product T=A'*Omega*A, where A is the arrow
282 : bidiagonal matrix formed by d, e. T is stored in dd, ee */
283 380 : PetscCall(DSAllocateWork_Private(ds,(n+6)*ld,4*ld,2*ld));
284 380 : R = ds->work+swu;
285 380 : swu += n*ld;
286 380 : perm = ds->iwork+iwu;
287 380 : iwu += n;
288 380 : cmplx = ds->iwork+iwu;
289 380 : dd = ds->rwork+rwu;
290 380 : rwu += ld;
291 380 : ee = ds->rwork+rwu;
292 380 : rwu += ld;
293 728 : for (i=0;i<l;i++) {dd[i] = d[i]*d[i]*Omega[i]; ee[i] = 0.0;}
294 5221 : for (i=l;i<=ds->k;i++) {
295 4841 : dd[i] = Omega[i]*d[i]*d[i];
296 4841 : ee[i] = Omega[i]*d[i]*e[i];
297 : }
298 4841 : for (i=l;i<k;i++) dd[k] += Omega[i]*e[i]*e[i];
299 4901 : for (i=k+1;i<n;i++) {
300 4521 : dd[i] = Omega[i]*d[i]*d[i]+Omega[i-1]*e[i-1]*e[i-1];
301 4521 : ee[i] = Omega[i]*d[i]*e[i];
302 : }
303 :
304 : /* Reduce T to tridiagonal form */
305 380 : PetscCall(DSArrowTridiag(n2,dd+l,ee+l,V+off,ld));
306 :
307 : /* Solve the tridiagonal eigenproblem corresponding to T */
308 380 : PetscCallBLAS("LAPACKsteqr",LAPACKsteqr_("V",&n1,dd+l,ee+l,V+off,&ld,ds->rwork+rwu,&info));
309 380 : SlepcCheckLapackInfo("steqr",info);
310 9742 : for (i=l;i<n;i++) wr[i] = PetscSqrtScalar(PetscAbs(dd[i]));
311 :
312 : /* Build left singular vectors: U=A*V*Sigma^-1 */
313 380 : PetscCall(PetscArrayzero(U+l*ld,n1*ld));
314 9362 : for (i=l;i<n-1;i++) {
315 8982 : scal = d[i];
316 8982 : PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+l*ld+i,&ld,U+l*ld+i,&ld));
317 8982 : j = (i<k)?k:i+1;
318 8982 : scal = e[i];
319 8982 : PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+l*ld+j,&ld,U+l*ld+i,&ld));
320 : }
321 380 : scal = d[n-1];
322 380 : PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+off+(n1-1),&ld,U+off+(n1-1),&ld));
323 : /* Multiply by Sigma^-1 */
324 9742 : for (i=l;i<n;i++) {scal = 1.0/wr[i]; PetscCallBLAS("BLASscal",BLASscal_(&n1,&scal,U+i*ld+l,&one));}
325 :
326 : } else { /* non-compact */
327 :
328 3 : PetscCall(DSAllocateWork_Private(ds,(n+6)*ld,PetscDefined(USE_COMPLEX)?4*ld:ld,2*ld));
329 3 : R = ds->work+swu;
330 3 : swu += n*ld;
331 3 : perm = ds->iwork+iwu;
332 3 : iwu += n;
333 3 : cmplx = ds->iwork+iwu;
334 3 : dd = ds->rwork+rwu;
335 533 : for (j=l;j<m;j++) {
336 371240 : for (i=0;i<n;i++) ds->work[i] = Omega[i]*A[i+j*ld];
337 530 : PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&m,&sone,A,&ld,ds->work,&incx,&szero,V+j*ld,&incx));
338 : }
339 :
340 : /* compute eigenvalues */
341 3 : lwork = (n+6)*ld;
342 : #if defined(PETSC_USE_COMPLEX)
343 : rwu += ld;
344 : PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&m,V,&ld,dd,ds->work,&lwork,ds->rwork+rwu,&info));
345 : #else
346 3 : PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&m,V,&ld,dd,ds->work,&lwork,&info));
347 : #endif
348 3 : SlepcCheckLapackInfo("syev",info);
349 533 : for (i=l;i<PetscMin(n,m);i++) d[i] = PetscSqrtReal(PetscAbsReal(dd[i]));
350 :
351 : /* Build left singular vectors: U=A*V*Sigma^-1 */
352 533 : for (j=l;j<PetscMin(n,m);j++) {
353 530 : scal = 1.0/d[j];
354 530 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&m,&scal,A,&ld,V+j*ld,&incx,&szero,U+j*ld,&incx));
355 : }
356 : }
357 :
358 383 : if (ctx->reorth) { /* Reinforce orthogonality */
359 2 : nv = n1;
360 22 : for (i=0;i<n;i++) cmplx[i] = 0;
361 2 : PetscCall(DSPseudoOrthog_HR(&nv,U+off,ld,Omega+l,R,ld,perm,cmplx,NULL,ds->work+swu));
362 : } else { /* Update Omega */
363 10257 : for (i=l;i<PetscMin(n,m);i++) Omega[i] = PetscSign(dd[i]);
364 : }
365 :
366 : /* Update projected problem */
367 383 : if (ds->compact) {
368 9742 : for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
369 380 : PetscCall(PetscArrayzero(e,n-1));
370 : } else {
371 533 : for (i=l;i<m;i++) PetscCall(PetscArrayzero(A+l+i*ld,n-l));
372 1251 : for (i=l;i<n;i++) A[i+i*ld] = d[i];
373 : }
374 10275 : for (i=l;i<PetscMin(n,m);i++) wr[i] = d[i];
375 383 : if (wi) for (i=l;i<PetscMin(n,m);i++) wi[i] = 0.0;
376 :
377 383 : if (ctx->reorth) { /* Update vectors V with R */
378 2 : scal = -1.0;
379 18 : for (i=0;i<nv;i++) {
380 16 : if (PetscRealPart(R[i+i*ld]) < 0.0) PetscCallBLAS("BLASscal",BLASscal_(&n1,&scal,V+(i+l)*ld+l,&one));
381 : }
382 : }
383 :
384 383 : if (!ds->compact) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
385 383 : PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_U],&U));
386 383 : PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_V],&V));
387 383 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
388 383 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&Omega));
389 383 : PetscFunctionReturn(PETSC_SUCCESS);
390 : }
391 :
392 : #if !defined(PETSC_HAVE_MPIUNI)
393 12 : static PetscErrorCode DSSynchronize_HSVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
394 : {
395 12 : PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
396 12 : PetscMPIInt n,rank,off=0,size,ldn,ld3,ld_;
397 12 : PetscScalar *A,*U,*V;
398 12 : PetscReal *T,*D;
399 :
400 12 : PetscFunctionBegin;
401 12 : if (ds->compact) kr = 3*ld;
402 0 : else k = (ds->n-l)*ld;
403 12 : kr += ld;
404 12 : if (ds->state>DS_STATE_RAW) k += 2*(ds->n-l)*ld;
405 12 : if (eigr) k += ds->n-l;
406 12 : PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
407 12 : PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
408 12 : PetscCall(PetscMPIIntCast(ds->n-l,&n));
409 12 : PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
410 12 : PetscCall(PetscMPIIntCast(3*ld,&ld3));
411 12 : PetscCall(PetscMPIIntCast(ld,&ld_));
412 12 : if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
413 0 : else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
414 12 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
415 12 : if (ds->state>DS_STATE_RAW) {
416 12 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
417 12 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
418 : }
419 12 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
420 12 : if (!rank) {
421 6 : if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
422 0 : else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
423 6 : PetscCallMPI(MPI_Pack(D,ld_,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
424 6 : if (ds->state>DS_STATE_RAW) {
425 6 : PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
426 6 : PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
427 : }
428 6 : if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
429 : }
430 24 : PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
431 12 : if (rank) {
432 6 : if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
433 0 : else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
434 6 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,D,ld_,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
435 6 : if (ds->state>DS_STATE_RAW) {
436 6 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
437 6 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
438 : }
439 6 : if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
440 : }
441 12 : if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
442 0 : else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
443 12 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
444 12 : if (ds->state>DS_STATE_RAW) {
445 12 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
446 12 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
447 : }
448 12 : PetscFunctionReturn(PETSC_SUCCESS);
449 : }
450 : #endif
451 :
452 1120 : static PetscErrorCode DSMatGetSize_HSVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
453 : {
454 1120 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
455 :
456 1120 : PetscFunctionBegin;
457 1120 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
458 1120 : switch (t) {
459 3 : case DS_MAT_A:
460 3 : *rows = ds->n;
461 3 : *cols = ds->extrarow? ctx->m+1: ctx->m;
462 3 : break;
463 0 : case DS_MAT_T:
464 0 : *rows = ds->n;
465 0 : *cols = PetscDefined(USE_COMPLEX)? 2: 3;
466 0 : break;
467 365 : case DS_MAT_D:
468 365 : *rows = ds->n;
469 365 : *cols = 1;
470 365 : break;
471 376 : case DS_MAT_U:
472 376 : *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
473 376 : *cols = ds->n;
474 376 : break;
475 376 : case DS_MAT_V:
476 376 : *rows = ds->state==DS_STATE_TRUNCATED? ctx->t: ctx->m;
477 376 : *cols = ctx->m;
478 376 : break;
479 0 : default:
480 0 : SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
481 : }
482 1120 : PetscFunctionReturn(PETSC_SUCCESS);
483 : }
484 :
485 383 : static PetscErrorCode DSHSVDSetDimensions_HSVD(DS ds,PetscInt m)
486 : {
487 383 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
488 :
489 383 : PetscFunctionBegin;
490 383 : DSCheckAlloc(ds,1);
491 383 : if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
492 0 : ctx->m = ds->ld;
493 : } else {
494 383 : PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
495 383 : ctx->m = m;
496 : }
497 383 : PetscFunctionReturn(PETSC_SUCCESS);
498 : }
499 :
500 : /*@
501 : DSHSVDSetDimensions - Sets the number of columns for a DSHSVD.
502 :
503 : Logically Collective
504 :
505 : Input Parameters:
506 : + ds - the direct solver context
507 : - m - the number of columns
508 :
509 : Notes:
510 : This call is complementary to DSSetDimensions(), to provide a dimension
511 : that is specific to this DS type.
512 :
513 : Level: intermediate
514 :
515 : .seealso: DSHSVDGetDimensions(), DSSetDimensions()
516 : @*/
517 383 : PetscErrorCode DSHSVDSetDimensions(DS ds,PetscInt m)
518 : {
519 383 : PetscFunctionBegin;
520 383 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
521 1149 : PetscValidLogicalCollectiveInt(ds,m,2);
522 383 : PetscTryMethod(ds,"DSHSVDSetDimensions_C",(DS,PetscInt),(ds,m));
523 383 : PetscFunctionReturn(PETSC_SUCCESS);
524 : }
525 :
526 1 : static PetscErrorCode DSHSVDGetDimensions_HSVD(DS ds,PetscInt *m)
527 : {
528 1 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
529 :
530 1 : PetscFunctionBegin;
531 1 : *m = ctx->m;
532 1 : PetscFunctionReturn(PETSC_SUCCESS);
533 : }
534 :
535 : /*@
536 : DSHSVDGetDimensions - Returns the number of columns for a DSHSVD.
537 :
538 : Not Collective
539 :
540 : Input Parameter:
541 : . ds - the direct solver context
542 :
543 : Output Parameters:
544 : . m - the number of columns
545 :
546 : Level: intermediate
547 :
548 : .seealso: DSHSVDSetDimensions()
549 : @*/
550 1 : PetscErrorCode DSHSVDGetDimensions(DS ds,PetscInt *m)
551 : {
552 1 : PetscFunctionBegin;
553 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
554 1 : PetscAssertPointer(m,2);
555 1 : PetscUseMethod(ds,"DSHSVDGetDimensions_C",(DS,PetscInt*),(ds,m));
556 1 : PetscFunctionReturn(PETSC_SUCCESS);
557 : }
558 :
559 4 : static PetscErrorCode DSHSVDSetReorthogonalize_HSVD(DS ds,PetscBool reorth)
560 : {
561 4 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
562 :
563 4 : PetscFunctionBegin;
564 4 : ctx->reorth = reorth;
565 4 : PetscFunctionReturn(PETSC_SUCCESS);
566 : }
567 :
568 : /*@
569 : DSHSVDSetReorthogonalize - Sets the reorthogonalization of the left vectors in a DSHSVD.
570 :
571 : Logically Collective
572 :
573 : Input Parameters:
574 : + ds - the direct solver context
575 : - reorth - the reorthogonalization flag
576 :
577 : Options Database Key:
578 : . -ds_hsvd_reorthog <bool> - sets the reorthogonalization flag
579 :
580 : Note:
581 : The computed left vectors (U) should be orthogonal with respect to the signature (D).
582 : But it may be necessary to enforce this with a final reorthogonalization step (omitted
583 : by default).
584 :
585 : Level: intermediate
586 :
587 : .seealso: DSHSVDGetReorthogonalize()
588 : @*/
589 4 : PetscErrorCode DSHSVDSetReorthogonalize(DS ds,PetscBool reorth)
590 : {
591 4 : PetscFunctionBegin;
592 4 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
593 12 : PetscValidLogicalCollectiveBool(ds,reorth,2);
594 4 : PetscTryMethod(ds,"DSHSVDSetReorthogonalize_C",(DS,PetscBool),(ds,reorth));
595 4 : PetscFunctionReturn(PETSC_SUCCESS);
596 : }
597 :
598 4 : static PetscErrorCode DSHSVDGetReorthogonalize_HSVD(DS ds,PetscBool *reorth)
599 : {
600 4 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
601 :
602 4 : PetscFunctionBegin;
603 4 : *reorth = ctx->reorth;
604 4 : PetscFunctionReturn(PETSC_SUCCESS);
605 : }
606 :
607 : /*@
608 : DSHSVDGetReorthogonalize - Returns the reorthogonalization flag of a DSHSVD.
609 :
610 : Not Collective
611 :
612 : Input Parameter:
613 : . ds - the direct solver context
614 :
615 : Output Parameters:
616 : . reorth - the reorthogonalization flag
617 :
618 : Level: intermediate
619 :
620 : .seealso: DSHSVDSetReorthogonalize()
621 : @*/
622 4 : PetscErrorCode DSHSVDGetReorthogonalize(DS ds,PetscBool *reorth)
623 : {
624 4 : PetscFunctionBegin;
625 4 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
626 4 : PetscAssertPointer(reorth,2);
627 4 : PetscUseMethod(ds,"DSHSVDGetReorthogonalize_C",(DS,PetscBool*),(ds,reorth));
628 4 : PetscFunctionReturn(PETSC_SUCCESS);
629 : }
630 :
631 20 : static PetscErrorCode DSSetFromOptions_HSVD(DS ds,PetscOptionItems *PetscOptionsObject)
632 : {
633 20 : PetscBool flg,reorth;
634 :
635 20 : PetscFunctionBegin;
636 20 : PetscOptionsHeadBegin(PetscOptionsObject,"DS HSVD Options");
637 :
638 20 : PetscCall(PetscOptionsBool("-ds_hsvd_reorthog","Reorthogonalize U vectors","DSHSVDSetReorthogonalize",PETSC_FALSE,&reorth,&flg));
639 20 : if (flg) PetscCall(DSHSVDSetReorthogonalize(ds,reorth));
640 :
641 20 : PetscOptionsHeadEnd();
642 20 : PetscFunctionReturn(PETSC_SUCCESS);
643 : }
644 :
645 21 : static PetscErrorCode DSDestroy_HSVD(DS ds)
646 : {
647 21 : PetscFunctionBegin;
648 21 : PetscCall(PetscFree(ds->data));
649 21 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetDimensions_C",NULL));
650 21 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetDimensions_C",NULL));
651 21 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetReorthogonalize_C",NULL));
652 21 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetReorthogonalize_C",NULL));
653 21 : PetscFunctionReturn(PETSC_SUCCESS);
654 : }
655 :
656 4 : static PetscErrorCode DSSetCompact_HSVD(DS ds,PetscBool comp)
657 : {
658 4 : PetscFunctionBegin;
659 4 : if (!comp) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
660 4 : PetscFunctionReturn(PETSC_SUCCESS);
661 : }
662 :
663 : /*MC
664 : DSHSVD - Dense Hyperbolic Singular Value Decomposition.
665 :
666 : Level: beginner
667 :
668 : Notes:
669 : The problem is expressed as A = U*Sigma*V', where A is rectangular in
670 : general, with n rows and m columns. U is orthogonal with respect to a
671 : signature matrix, stored in D. V is orthogonal. Sigma is a diagonal
672 : matrix whose diagonal elements are the arguments of DSSolve(). After
673 : solve, A is overwritten with Sigma, D is overwritten with the new signature.
674 :
675 : The matrices of left and right singular vectors, U and V, have size n and m,
676 : respectively. The number of columns m must be specified via DSHSVDSetDimensions().
677 :
678 : If the DS object is in the intermediate state, A is assumed to be in upper
679 : bidiagonal form (possibly with an arrow) and is stored in compact format
680 : on matrix T. The compact storage is implemented for the square case
681 : only, m=n. The extra row should be interpreted in this case as an extra column.
682 :
683 : Used DS matrices:
684 : + DS_MAT_A - problem matrix (used only if compact=false)
685 : . DS_MAT_T - upper bidiagonal matrix
686 : . DS_MAT_D - diagonal matrix (signature)
687 : . DS_MAT_U - left singular vectors
688 : - DS_MAT_V - right singular vectors
689 :
690 : Implemented methods:
691 : . 0 - Cross product A'*Omega*A
692 :
693 : .seealso: DSCreate(), DSSetType(), DSType, DSHSVDSetDimensions()
694 : M*/
695 21 : SLEPC_EXTERN PetscErrorCode DSCreate_HSVD(DS ds)
696 : {
697 21 : DS_HSVD *ctx;
698 :
699 21 : PetscFunctionBegin;
700 21 : PetscCall(PetscNew(&ctx));
701 21 : ds->data = (void*)ctx;
702 :
703 21 : ds->ops->allocate = DSAllocate_HSVD;
704 21 : ds->ops->setfromoptions = DSSetFromOptions_HSVD;
705 21 : ds->ops->view = DSView_HSVD;
706 21 : ds->ops->vectors = DSVectors_HSVD;
707 21 : ds->ops->solve[0] = DSSolve_HSVD_CROSS;
708 21 : ds->ops->sort = DSSort_HSVD;
709 21 : ds->ops->truncate = DSTruncate_HSVD;
710 21 : ds->ops->update = DSUpdateExtraRow_HSVD;
711 21 : ds->ops->destroy = DSDestroy_HSVD;
712 21 : ds->ops->matgetsize = DSMatGetSize_HSVD;
713 : #if !defined(PETSC_HAVE_MPIUNI)
714 21 : ds->ops->synchronize = DSSynchronize_HSVD;
715 : #endif
716 21 : ds->ops->setcompact = DSSetCompact_HSVD;
717 21 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetDimensions_C",DSHSVDSetDimensions_HSVD));
718 21 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetDimensions_C",DSHSVDGetDimensions_HSVD));
719 21 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetReorthogonalize_C",DSHSVDSetReorthogonalize_HSVD));
720 21 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetReorthogonalize_C",DSHSVDGetReorthogonalize_HSVD));
721 21 : PetscFunctionReturn(PETSC_SUCCESS);
722 : }
|
