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 18 : static PetscErrorCode DSAllocate_HSVD(DS ds,PetscInt ld)
21 : {
22 18 : PetscFunctionBegin;
23 18 : if (!ds->compact) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
24 18 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
25 18 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_V));
26 18 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
27 18 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_D));
28 18 : PetscCall(PetscFree(ds->perm));
29 18 : PetscCall(PetscMalloc1(ld,&ds->perm));
30 18 : 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 97 : static PetscErrorCode DSSort_HSVD(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
154 : {
155 97 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
156 97 : PetscInt n,l,i,*perm,ld=ds->ld;
157 97 : PetscScalar *A;
158 97 : PetscReal *d,*s;
159 :
160 97 : PetscFunctionBegin;
161 97 : if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
162 97 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
163 97 : l = ds->l;
164 97 : n = PetscMin(ds->n,ctx->m);
165 97 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
166 97 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
167 97 : PetscCall(DSAllocateWork_Private(ds,0,ds->ld,0));
168 97 : perm = ds->perm;
169 97 : if (!rr) PetscCall(DSSortEigenvaluesReal_Private(ds,d,perm));
170 0 : else PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
171 97 : PetscCall(PetscArraycpy(ds->rwork,s,n));
172 2025 : for (i=l;i<n;i++) s[i] = ds->rwork[perm[i]];
173 2025 : for (i=l;i<n;i++) wr[i] = d[perm[i]];
174 97 : PetscCall(DSPermuteBoth_Private(ds,l,n,ds->n,ctx->m,DS_MAT_U,DS_MAT_V,perm));
175 2025 : for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
176 97 : if (!ds->compact) {
177 2 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
178 332 : for (i=l;i<n;i++) A[i+i*ld] = wr[i];
179 2 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
180 : }
181 97 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
182 97 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
183 97 : PetscFunctionReturn(PETSC_SUCCESS);
184 : }
185 :
186 93 : static PetscErrorCode DSUpdateExtraRow_HSVD(DS ds)
187 : {
188 93 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
189 93 : PetscInt i;
190 93 : PetscBLASInt n=0,m=0,ld,l;
191 93 : const PetscScalar *U;
192 93 : PetscReal *T,*e,*Omega,beta;
193 :
194 93 : PetscFunctionBegin;
195 93 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
196 93 : PetscCall(PetscBLASIntCast(ds->n,&n));
197 93 : PetscCall(PetscBLASIntCast(ctx->m,&m));
198 93 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
199 93 : PetscCall(PetscBLASIntCast(ds->l,&l));
200 93 : PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_U],&U));
201 93 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&Omega));
202 93 : PetscCheck(ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-compact storage");
203 93 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
204 93 : e = T+ld;
205 93 : beta = PetscAbs(e[m-1]); /* in compact, we assume all entries are zero except the last one */
206 2001 : for (i=0;i<n;i++) e[i] = PetscRealPart(beta*U[n-1+i*ld]*Omega[i]);
207 93 : ds->k = m;
208 93 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
209 93 : PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_U],&U));
210 93 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&Omega));
211 93 : PetscFunctionReturn(PETSC_SUCCESS);
212 : }
213 :
214 91 : static PetscErrorCode DSTruncate_HSVD(DS ds,PetscInt n,PetscBool trim)
215 : {
216 91 : PetscInt i,ld=ds->ld,l=ds->l;
217 91 : PetscScalar *A;
218 91 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
219 :
220 91 : PetscFunctionBegin;
221 91 : if (!ds->compact && ds->extrarow) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
222 91 : if (trim) {
223 12 : 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 12 : ds->l = 0;
227 12 : ds->k = 0;
228 12 : ds->n = n;
229 12 : ctx->m = n;
230 12 : ds->t = ds->n; /* truncated length equal to the new dimension */
231 12 : ctx->t = ctx->m; /* must also keep the previous dimension of V */
232 : } else {
233 79 : 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 79 : ds->k = (ds->extrarow)? n: 0;
239 79 : ds->t = ds->n; /* truncated length equal to previous dimension */
240 79 : ctx->t = ctx->m; /* must also keep the previous dimension of V */
241 79 : ds->n = n;
242 79 : ctx->m = n;
243 : }
244 91 : if (!ds->compact && ds->extrarow) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
245 91 : PetscFunctionReturn(PETSC_SUCCESS);
246 : }
247 :
248 97 : static PetscErrorCode DSSolve_HSVD_CROSS(DS ds,PetscScalar *wr,PetscScalar *wi)
249 : {
250 97 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
251 97 : PetscInt i,j,k=ds->k,rwu=0,iwu=0,swu=0,nv;
252 97 : PetscBLASInt n1,n2,info,l=0,n=0,m=0,ld,off,one=1,*perm,*cmplx,incx=1,lwork;
253 97 : PetscScalar *A,*U,*V,scal,*R,sone=1.0,szero=0.0;
254 97 : PetscReal *d,*e,*dd,*ee,*Omega;
255 :
256 97 : PetscFunctionBegin;
257 97 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
258 97 : PetscCall(PetscBLASIntCast(ds->n,&n));
259 97 : PetscCall(PetscBLASIntCast(ctx->m,&m));
260 97 : PetscCheck(!ds->compact || n==m,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-square matrices in compact storage");
261 97 : PetscCheck(ds->compact || n>=m,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for the case of more columns than rows");
262 97 : PetscCall(PetscBLASIntCast(ds->l,&l));
263 97 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
264 97 : PetscCall(PetscBLASIntCast(PetscMax(0,ds->k-ds->l+1),&n2));
265 97 : n1 = n-l; /* n1 = size of leading block, excl. locked + size of trailing block */
266 97 : off = l+l*ld;
267 97 : if (!ds->compact) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
268 97 : PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_U],&U));
269 97 : PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_V],&V));
270 97 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
271 97 : e = d+ld;
272 97 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&Omega));
273 97 : PetscCall(PetscArrayzero(U,ld*ld));
274 427 : for (i=0;i<l;i++) U[i+i*ld] = 1.0;
275 97 : PetscCall(PetscArrayzero(V,ld*ld));
276 3073 : for (i=0;i<n;i++) V[i+i*ld] = 1.0;
277 427 : for (i=0;i<l;i++) wr[i] = d[i];
278 97 : if (wi) for (i=0;i<l;i++) wi[i] = 0.0;
279 :
280 97 : 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 95 : PetscCall(DSAllocateWork_Private(ds,(n+6)*ld,4*ld,2*ld));
284 95 : R = ds->work+swu;
285 95 : swu += n*ld;
286 95 : perm = ds->iwork+iwu;
287 95 : iwu += n;
288 95 : cmplx = ds->iwork+iwu;
289 95 : dd = ds->rwork+rwu;
290 95 : rwu += ld;
291 95 : ee = ds->rwork+rwu;
292 95 : rwu += ld;
293 425 : for (i=0;i<l;i++) {dd[i] = d[i]*d[i]*Omega[i]; ee[i] = 0.0;}
294 789 : for (i=l;i<=ds->k;i++) {
295 694 : dd[i] = Omega[i]*d[i]*d[i];
296 694 : ee[i] = Omega[i]*d[i]*e[i];
297 : }
298 694 : for (i=l;i<k;i++) dd[k] += Omega[i]*e[i]*e[i];
299 999 : for (i=k+1;i<n;i++) {
300 904 : dd[i] = Omega[i]*d[i]*d[i]+Omega[i-1]*e[i-1]*e[i-1];
301 904 : ee[i] = Omega[i]*d[i]*e[i];
302 : }
303 :
304 : /* Reduce T to tridiagonal form */
305 95 : PetscCall(DSArrowTridiag(n2,dd+l,ee+l,V+off,ld));
306 :
307 : /* Solve the tridiagonal eigenproblem corresponding to T */
308 95 : PetscCallBLAS("LAPACKsteqr",LAPACKsteqr_("V",&n1,dd+l,ee+l,V+off,&ld,ds->rwork+rwu,&info));
309 95 : SlepcCheckLapackInfo("steqr",info);
310 1693 : for (i=l;i<n;i++) wr[i] = PetscSqrtScalar(PetscAbs(dd[i]));
311 :
312 : /* Build left singular vectors: U=A*V*Sigma^-1 */
313 95 : PetscCall(PetscArrayzero(U+l*ld,n1*ld));
314 1598 : for (i=l;i<n-1;i++) {
315 1503 : scal = d[i];
316 1503 : PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+l*ld+i,&ld,U+l*ld+i,&ld));
317 1503 : j = (i<k)?k:i+1;
318 1503 : scal = e[i];
319 1503 : PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+l*ld+j,&ld,U+l*ld+i,&ld));
320 : }
321 95 : scal = d[n-1];
322 95 : PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+off+(n1-1),&ld,U+off+(n1-1),&ld));
323 : /* Multiply by Sigma^-1 */
324 1693 : 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 2 : PetscCall(DSAllocateWork_Private(ds,(n+6)*ld,PetscDefined(USE_COMPLEX)?4*ld:ld,2*ld));
329 2 : R = ds->work+swu;
330 2 : swu += n*ld;
331 2 : perm = ds->iwork+iwu;
332 2 : iwu += n;
333 2 : cmplx = ds->iwork+iwu;
334 2 : dd = ds->rwork+rwu;
335 332 : for (j=l;j<m;j++) {
336 331040 : for (i=0;i<n;i++) ds->work[i] = Omega[i]*A[i+j*ld];
337 330 : PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&m,&sone,A,&ld,ds->work,&incx,&szero,V+j*ld,&incx));
338 : }
339 :
340 : /* compute eigenvalues */
341 2 : lwork = (n+6)*ld;
342 : #if defined(PETSC_USE_COMPLEX)
343 2 : rwu += ld;
344 2 : PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&m,V,&ld,dd,ds->work,&lwork,ds->rwork+rwu,&info));
345 : #else
346 : PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&m,V,&ld,dd,ds->work,&lwork,&info));
347 : #endif
348 2 : SlepcCheckLapackInfo("syev",info);
349 332 : 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 332 : for (j=l;j<PetscMin(n,m);j++) {
353 330 : scal = 1.0/d[j];
354 330 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&m,&scal,A,&ld,V+j*ld,&incx,&szero,U+j*ld,&incx));
355 : }
356 : }
357 :
358 97 : 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 2007 : for (i=l;i<PetscMin(n,m);i++) Omega[i] = PetscSign(dd[i]);
364 : }
365 :
366 : /* Update projected problem */
367 97 : if (ds->compact) {
368 1693 : for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
369 95 : PetscCall(PetscArrayzero(e,n-1));
370 : } else {
371 332 : for (i=l;i<m;i++) PetscCall(PetscArrayzero(A+l+i*ld,n-l));
372 1050 : for (i=l;i<n;i++) A[i+i*ld] = d[i];
373 : }
374 2025 : for (i=l;i<PetscMin(n,m);i++) wr[i] = d[i];
375 97 : if (wi) for (i=l;i<PetscMin(n,m);i++) wi[i] = 0.0;
376 :
377 97 : 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 97 : if (!ds->compact) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
385 97 : PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_U],&U));
386 97 : PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_V],&V));
387 97 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
388 97 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&Omega));
389 97 : PetscFunctionReturn(PETSC_SUCCESS);
390 : }
391 :
392 : #if !defined(PETSC_HAVE_MPIUNI)
393 14 : static PetscErrorCode DSSynchronize_HSVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
394 : {
395 14 : PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
396 14 : PetscMPIInt n,rank,off=0,size,ldn,ld3,ld_;
397 14 : PetscScalar *A,*U,*V;
398 14 : PetscReal *T,*D;
399 :
400 14 : PetscFunctionBegin;
401 14 : if (ds->compact) kr = 3*ld;
402 0 : else k = (ds->n-l)*ld;
403 14 : kr += ld;
404 14 : if (ds->state>DS_STATE_RAW) k += 2*(ds->n-l)*ld;
405 14 : if (eigr) k += ds->n-l;
406 14 : PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
407 14 : PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
408 14 : PetscCall(PetscMPIIntCast(ds->n-l,&n));
409 14 : PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
410 14 : PetscCall(PetscMPIIntCast(3*ld,&ld3));
411 14 : PetscCall(PetscMPIIntCast(ld,&ld_));
412 14 : if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
413 0 : else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
414 14 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
415 14 : if (ds->state>DS_STATE_RAW) {
416 14 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
417 14 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
418 : }
419 14 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
420 14 : if (!rank) {
421 7 : 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 7 : PetscCallMPI(MPI_Pack(D,ld_,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
424 7 : if (ds->state>DS_STATE_RAW) {
425 7 : PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
426 7 : PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
427 : }
428 7 : if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
429 : }
430 28 : PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
431 14 : if (rank) {
432 7 : 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 7 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,D,ld_,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
435 7 : if (ds->state>DS_STATE_RAW) {
436 7 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
437 7 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
438 : }
439 7 : if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
440 : }
441 14 : if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
442 0 : else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
443 14 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
444 14 : if (ds->state>DS_STATE_RAW) {
445 14 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
446 14 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
447 : }
448 14 : PetscFunctionReturn(PETSC_SUCCESS);
449 : }
450 : #endif
451 :
452 265 : static PetscErrorCode DSMatGetSize_HSVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
453 : {
454 265 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
455 :
456 265 : PetscFunctionBegin;
457 265 : PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
458 265 : switch (t) {
459 2 : case DS_MAT_A:
460 2 : *rows = ds->n;
461 2 : *cols = ds->extrarow? ctx->m+1: ctx->m;
462 2 : break;
463 0 : case DS_MAT_T:
464 0 : *rows = ds->n;
465 0 : *cols = PetscDefined(USE_COMPLEX)? 2: 3;
466 0 : break;
467 81 : case DS_MAT_D:
468 81 : *rows = ds->n;
469 81 : *cols = 1;
470 81 : break;
471 91 : case DS_MAT_U:
472 91 : *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
473 91 : *cols = ds->n;
474 91 : break;
475 91 : case DS_MAT_V:
476 91 : *rows = ds->state==DS_STATE_TRUNCATED? ctx->t: ctx->m;
477 91 : *cols = ctx->m;
478 91 : break;
479 0 : default:
480 0 : SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
481 : }
482 265 : PetscFunctionReturn(PETSC_SUCCESS);
483 : }
484 :
485 97 : static PetscErrorCode DSHSVDSetDimensions_HSVD(DS ds,PetscInt m)
486 : {
487 97 : DS_HSVD *ctx = (DS_HSVD*)ds->data;
488 :
489 97 : PetscFunctionBegin;
490 97 : DSCheckAlloc(ds,1);
491 97 : if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
492 0 : ctx->m = ds->ld;
493 : } else {
494 97 : PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
495 97 : ctx->m = m;
496 : }
497 97 : 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 97 : PetscErrorCode DSHSVDSetDimensions(DS ds,PetscInt m)
518 : {
519 97 : PetscFunctionBegin;
520 97 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
521 291 : PetscValidLogicalCollectiveInt(ds,m,2);
522 97 : PetscTryMethod(ds,"DSHSVDSetDimensions_C",(DS,PetscInt),(ds,m));
523 97 : 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 17 : static PetscErrorCode DSSetFromOptions_HSVD(DS ds,PetscOptionItems *PetscOptionsObject)
632 : {
633 17 : PetscBool flg,reorth;
634 :
635 17 : PetscFunctionBegin;
636 17 : PetscOptionsHeadBegin(PetscOptionsObject,"DS HSVD Options");
637 :
638 17 : PetscCall(PetscOptionsBool("-ds_hsvd_reorthog","Reorthogonalize U vectors","DSHSVDSetReorthogonalize",PETSC_FALSE,&reorth,&flg));
639 17 : if (flg) PetscCall(DSHSVDSetReorthogonalize(ds,reorth));
640 :
641 17 : PetscOptionsHeadEnd();
642 17 : PetscFunctionReturn(PETSC_SUCCESS);
643 : }
644 :
645 18 : static PetscErrorCode DSDestroy_HSVD(DS ds)
646 : {
647 18 : PetscFunctionBegin;
648 18 : PetscCall(PetscFree(ds->data));
649 18 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetDimensions_C",NULL));
650 18 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetDimensions_C",NULL));
651 18 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetReorthogonalize_C",NULL));
652 18 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetReorthogonalize_C",NULL));
653 18 : 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 1 : static PetscErrorCode DSReallocate_HSVD(DS ds,PetscInt ld)
664 : {
665 1 : PetscInt i,*perm=ds->perm;
666 :
667 1 : PetscFunctionBegin;
668 23 : for (i=0;i<DS_NUM_MAT;i++) {
669 22 : if (!ds->compact && i==DS_MAT_A) continue;
670 22 : if (i!=DS_MAT_U && i!=DS_MAT_V && i!=DS_MAT_T && i!=DS_MAT_D) PetscCall(MatDestroy(&ds->omat[i]));
671 : }
672 :
673 1 : if (!ds->compact) PetscCall(DSReallocateMat_Private(ds,DS_MAT_A,ld));
674 1 : PetscCall(DSReallocateMat_Private(ds,DS_MAT_U,ld));
675 1 : PetscCall(DSReallocateMat_Private(ds,DS_MAT_V,ld));
676 1 : PetscCall(DSReallocateMat_Private(ds,DS_MAT_T,ld));
677 1 : PetscCall(DSReallocateMat_Private(ds,DS_MAT_D,ld));
678 :
679 1 : PetscCall(PetscMalloc1(ld,&ds->perm));
680 1 : PetscCall(PetscArraycpy(ds->perm,perm,ds->ld));
681 1 : PetscCall(PetscFree(perm));
682 1 : PetscFunctionReturn(PETSC_SUCCESS);
683 : }
684 :
685 : /*MC
686 : DSHSVD - Dense Hyperbolic Singular Value Decomposition.
687 :
688 : Level: beginner
689 :
690 : Notes:
691 : The problem is expressed as A = U*Sigma*V', where A is rectangular in
692 : general, with n rows and m columns. U is orthogonal with respect to a
693 : signature matrix, stored in D. V is orthogonal. Sigma is a diagonal
694 : matrix whose diagonal elements are the arguments of DSSolve(). After
695 : solve, A is overwritten with Sigma, D is overwritten with the new signature.
696 :
697 : The matrices of left and right singular vectors, U and V, have size n and m,
698 : respectively. The number of columns m must be specified via DSHSVDSetDimensions().
699 :
700 : If the DS object is in the intermediate state, A is assumed to be in upper
701 : bidiagonal form (possibly with an arrow) and is stored in compact format
702 : on matrix T. The compact storage is implemented for the square case
703 : only, m=n. The extra row should be interpreted in this case as an extra column.
704 :
705 : Used DS matrices:
706 : + DS_MAT_A - problem matrix (used only if compact=false)
707 : . DS_MAT_T - upper bidiagonal matrix
708 : . DS_MAT_D - diagonal matrix (signature)
709 : . DS_MAT_U - left singular vectors
710 : - DS_MAT_V - right singular vectors
711 :
712 : Implemented methods:
713 : . 0 - Cross product A'*Omega*A
714 :
715 : .seealso: DSCreate(), DSSetType(), DSType, DSHSVDSetDimensions()
716 : M*/
717 18 : SLEPC_EXTERN PetscErrorCode DSCreate_HSVD(DS ds)
718 : {
719 18 : DS_HSVD *ctx;
720 :
721 18 : PetscFunctionBegin;
722 18 : PetscCall(PetscNew(&ctx));
723 18 : ds->data = (void*)ctx;
724 :
725 18 : ds->ops->allocate = DSAllocate_HSVD;
726 18 : ds->ops->setfromoptions = DSSetFromOptions_HSVD;
727 18 : ds->ops->view = DSView_HSVD;
728 18 : ds->ops->vectors = DSVectors_HSVD;
729 18 : ds->ops->solve[0] = DSSolve_HSVD_CROSS;
730 18 : ds->ops->sort = DSSort_HSVD;
731 18 : ds->ops->truncate = DSTruncate_HSVD;
732 18 : ds->ops->update = DSUpdateExtraRow_HSVD;
733 18 : ds->ops->destroy = DSDestroy_HSVD;
734 18 : ds->ops->matgetsize = DSMatGetSize_HSVD;
735 : #if !defined(PETSC_HAVE_MPIUNI)
736 18 : ds->ops->synchronize = DSSynchronize_HSVD;
737 : #endif
738 18 : ds->ops->setcompact = DSSetCompact_HSVD;
739 18 : ds->ops->reallocate = DSReallocate_HSVD;
740 18 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetDimensions_C",DSHSVDSetDimensions_HSVD));
741 18 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetDimensions_C",DSHSVDGetDimensions_HSVD));
742 18 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetReorthogonalize_C",DSHSVDSetReorthogonalize_HSVD));
743 18 : PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetReorthogonalize_C",DSHSVDGetReorthogonalize_HSVD));
744 18 : PetscFunctionReturn(PETSC_SUCCESS);
745 : }
|
