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