Actual source code: dssvd.c
slepc-3.21.1 2024-04-26
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: #include <slepc/private/dsimpl.h>
12: #include <slepcblaslapack.h>
14: typedef struct {
15: PetscInt m; /* number of columns */
16: PetscInt t; /* number of rows of V after truncating */
17: } DS_SVD;
19: static PetscErrorCode DSAllocate_SVD(DS ds,PetscInt ld)
20: {
21: PetscFunctionBegin;
22: PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
23: PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
24: PetscCall(DSAllocateMat_Private(ds,DS_MAT_V));
25: PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
26: PetscCall(PetscFree(ds->perm));
27: PetscCall(PetscMalloc1(ld,&ds->perm));
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
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: */
56: static PetscErrorCode DSSwitchFormat_SVD(DS ds)
57: {
58: DS_SVD *ctx = (DS_SVD*)ds->data;
59: PetscReal *T;
60: PetscScalar *A;
61: PetscInt i,m=ctx->m,k=ds->k,ld=ds->ld;
63: PetscFunctionBegin;
64: PetscCheck(m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
65: /* switch from compact (arrow) to dense storage */
66: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_A],&A));
67: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
68: PetscCall(PetscArrayzero(A,ld*ld));
69: for (i=0;i<k;i++) {
70: A[i+i*ld] = T[i];
71: A[i+k*ld] = T[i+ld];
72: }
73: A[k+k*ld] = T[k];
74: for (i=k+1;i<m;i++) {
75: A[i+i*ld] = T[i];
76: A[i-1+i*ld] = T[i-1+ld];
77: }
78: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_A],&A));
79: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: static PetscErrorCode DSView_SVD(DS ds,PetscViewer viewer)
84: {
85: DS_SVD *ctx = (DS_SVD*)ds->data;
86: PetscViewerFormat format;
87: PetscInt i,j,r,c,m=ctx->m,rows,cols;
88: PetscReal *T,value;
90: PetscFunctionBegin;
91: PetscCall(PetscViewerGetFormat(viewer,&format));
92: if (format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
93: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
94: PetscCall(PetscViewerASCIIPrintf(viewer,"number of columns: %" PetscInt_FMT "\n",m));
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
97: PetscCheck(m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
98: if (ds->compact) {
99: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
100: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
101: rows = ds->n;
102: cols = ds->extrarow? m+1: m;
103: if (format == PETSC_VIEWER_ASCII_MATLAB) {
104: PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rows,cols));
105: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",2*ds->n));
106: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
107: 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]));
108: for (i=0;i<cols-1;i++) {
109: r = PetscMax(i+2,ds->k+1);
110: c = i+1;
111: PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",c,r,(double)T[i+ds->ld]));
112: }
113: PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_T]));
114: } else {
115: for (i=0;i<rows;i++) {
116: for (j=0;j<cols;j++) {
117: if (i==j) value = T[i];
118: else if (i<ds->k && j==ds->k) value = T[PetscMin(i,j)+ds->ld];
119: else if (i+1==j && i>=ds->k) value = T[i+ds->ld];
120: else value = 0.0;
121: PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
122: }
123: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
124: }
125: }
126: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
127: PetscCall(PetscViewerFlush(viewer));
128: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
129: } else PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
130: if (ds->state>DS_STATE_INTERMEDIATE) {
131: PetscCall(DSViewMat(ds,viewer,DS_MAT_U));
132: PetscCall(DSViewMat(ds,viewer,DS_MAT_V));
133: }
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: static PetscErrorCode DSVectors_SVD(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
138: {
139: PetscFunctionBegin;
140: switch (mat) {
141: case DS_MAT_U:
142: case DS_MAT_V:
143: if (rnorm) *rnorm = 0.0;
144: break;
145: default:
146: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
147: }
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: static PetscErrorCode DSSort_SVD(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
152: {
153: DS_SVD *ctx = (DS_SVD*)ds->data;
154: PetscInt n,l,i,*perm,ld=ds->ld;
155: PetscScalar *A;
156: PetscReal *d;
158: PetscFunctionBegin;
159: if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
160: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
161: l = ds->l;
162: n = PetscMin(ds->n,ctx->m);
163: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
164: perm = ds->perm;
165: if (!rr) PetscCall(DSSortEigenvaluesReal_Private(ds,d,perm));
166: else PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
167: for (i=l;i<n;i++) wr[i] = d[perm[i]];
168: PetscCall(DSPermuteBoth_Private(ds,l,n,ds->n,ctx->m,DS_MAT_U,DS_MAT_V,perm));
169: for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
170: if (!ds->compact) {
171: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
172: for (i=l;i<n;i++) A[i+i*ld] = wr[i];
173: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
174: }
175: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: static PetscErrorCode DSUpdateExtraRow_SVD(DS ds)
180: {
181: DS_SVD *ctx = (DS_SVD*)ds->data;
182: PetscInt i;
183: PetscBLASInt n=0,m=0,ld,incx=1;
184: PetscScalar *A,*x,*y,one=1.0,zero=0.0;
185: PetscReal *T,*e,beta;
186: const PetscScalar *U;
188: PetscFunctionBegin;
189: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
190: PetscCall(PetscBLASIntCast(ds->n,&n));
191: PetscCall(PetscBLASIntCast(ctx->m,&m));
192: PetscCall(PetscBLASIntCast(ds->ld,&ld));
193: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_U],&U));
194: if (ds->compact) {
195: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
196: e = T+ld;
197: beta = e[m-1]; /* in compact, we assume all entries are zero except the last one */
198: for (i=0;i<n;i++) e[i] = PetscRealPart(beta*U[n-1+i*ld]);
199: ds->k = m;
200: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
201: } else {
202: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
203: PetscCall(DSAllocateWork_Private(ds,2*ld,0,0));
204: x = ds->work;
205: y = ds->work+ld;
206: for (i=0;i<n;i++) x[i] = PetscConj(A[i+m*ld]);
207: PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,U,&ld,x,&incx,&zero,y,&incx));
208: for (i=0;i<n;i++) A[i+m*ld] = PetscConj(y[i]);
209: ds->k = m;
210: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
211: }
212: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_U],&U));
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: static PetscErrorCode DSTruncate_SVD(DS ds,PetscInt n,PetscBool trim)
217: {
218: PetscInt i,ld=ds->ld,l=ds->l;
219: PetscScalar *A;
220: DS_SVD *ctx = (DS_SVD*)ds->data;
222: PetscFunctionBegin;
223: if (!ds->compact && ds->extrarow) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
224: if (trim) {
225: if (!ds->compact && ds->extrarow) { /* clean extra column */
226: for (i=l;i<ds->n;i++) A[i+ctx->m*ld] = 0.0;
227: }
228: ds->l = 0;
229: ds->k = 0;
230: ds->n = n;
231: ctx->m = n;
232: ds->t = ds->n; /* truncated length equal to the new dimension */
233: ctx->t = ctx->m; /* must also keep the previous dimension of V */
234: } else {
235: if (!ds->compact && ds->extrarow && ds->k==ds->n) {
236: /* copy entries of extra column to the new position, then clean last row */
237: for (i=l;i<n;i++) A[i+n*ld] = A[i+ctx->m*ld];
238: for (i=l;i<ds->n;i++) A[i+ctx->m*ld] = 0.0;
239: }
240: ds->k = (ds->extrarow)? n: 0;
241: ds->t = ds->n; /* truncated length equal to previous dimension */
242: ctx->t = ctx->m; /* must also keep the previous dimension of V */
243: ds->n = n;
244: ctx->m = n;
245: }
246: if (!ds->compact && ds->extrarow) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: static PetscErrorCode DSSolve_SVD_DC(DS ds,PetscScalar *wr,PetscScalar *wi)
251: {
252: DS_SVD *ctx = (DS_SVD*)ds->data;
253: PetscInt i,j;
254: PetscBLASInt n1,m1,info,l = 0,n = 0,m = 0,nm,ld,off,lwork;
255: PetscScalar *A,*U,*V,*W,qwork;
256: PetscReal *d,*e,*Ur,*Vr;
258: PetscFunctionBegin;
259: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
260: PetscCall(PetscBLASIntCast(ds->n,&n));
261: PetscCall(PetscBLASIntCast(ctx->m,&m));
262: PetscCall(PetscBLASIntCast(ds->l,&l));
263: PetscCall(PetscBLASIntCast(ds->ld,&ld));
264: n1 = n-l; /* n1 = size of leading block, excl. locked + size of trailing block */
265: m1 = m-l;
266: off = l+l*ld;
267: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
268: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_U],&U));
269: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_V],&V));
270: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
271: e = d+ld;
272: PetscCall(PetscArrayzero(U,ld*ld));
273: for (i=0;i<l;i++) U[i+i*ld] = 1.0;
274: PetscCall(PetscArrayzero(V,ld*ld));
275: for (i=0;i<l;i++) V[i+i*ld] = 1.0;
277: if (ds->state>DS_STATE_RAW) {
278: /* solve bidiagonal SVD problem */
279: for (i=0;i<l;i++) wr[i] = d[i];
280: #if defined(PETSC_USE_COMPLEX)
281: PetscCall(DSAllocateWork_Private(ds,0,3*n1*n1+4*n1+2*ld*ld,8*n1));
282: Ur = ds->rwork+3*n1*n1+4*n1;
283: Vr = ds->rwork+3*n1*n1+4*n1+ld*ld;
284: #else
285: PetscCall(DSAllocateWork_Private(ds,0,3*n1*n1+4*n1+ld*ld,8*n1));
286: Ur = U;
287: Vr = ds->rwork+3*n1*n1+4*n1;
288: #endif
289: PetscCallBLAS("LAPACKbdsdc",LAPACKbdsdc_("U","I",&n1,d+l,e+l,Ur+off,&ld,Vr+off,&ld,NULL,NULL,ds->rwork,ds->iwork,&info));
290: SlepcCheckLapackInfo("bdsdc",info);
291: for (i=l;i<n;i++) {
292: for (j=l;j<n;j++) {
293: #if defined(PETSC_USE_COMPLEX)
294: U[i+j*ld] = Ur[i+j*ld];
295: #endif
296: V[i+j*ld] = PetscConj(Vr[j+i*ld]); /* transpose VT returned by Lapack */
297: }
298: }
299: } else {
300: /* solve general rectangular SVD problem */
301: PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
302: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&W));
303: if (ds->compact) PetscCall(DSSwitchFormat_SVD(ds));
304: for (i=0;i<l;i++) wr[i] = d[i];
305: nm = PetscMin(n,m);
306: PetscCall(DSAllocateWork_Private(ds,0,0,8*nm));
307: lwork = -1;
308: #if defined(PETSC_USE_COMPLEX)
309: PetscCall(DSAllocateWork_Private(ds,0,5*nm*nm+7*nm,0));
310: PetscCallBLAS("LAPACKgesdd",LAPACKgesdd_("A",&n1,&m1,A+off,&ld,d+l,U+off,&ld,W+off,&ld,&qwork,&lwork,ds->rwork,ds->iwork,&info));
311: #else
312: PetscCallBLAS("LAPACKgesdd",LAPACKgesdd_("A",&n1,&m1,A+off,&ld,d+l,U+off,&ld,W+off,&ld,&qwork,&lwork,ds->iwork,&info));
313: #endif
314: SlepcCheckLapackInfo("gesdd",info);
315: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(qwork),&lwork));
316: PetscCall(DSAllocateWork_Private(ds,lwork,0,0));
317: #if defined(PETSC_USE_COMPLEX)
318: 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));
319: #else
320: PetscCallBLAS("LAPACKgesdd",LAPACKgesdd_("A",&n1,&m1,A+off,&ld,d+l,U+off,&ld,W+off,&ld,ds->work,&lwork,ds->iwork,&info));
321: #endif
322: SlepcCheckLapackInfo("gesdd",info);
323: for (i=l;i<m;i++) {
324: for (j=l;j<m;j++) V[i+j*ld] = PetscConj(W[j+i*ld]); /* transpose VT returned by Lapack */
325: }
326: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&W));
327: }
328: for (i=l;i<PetscMin(ds->n,ctx->m);i++) wr[i] = d[i];
330: /* create diagonal matrix as a result */
331: if (ds->compact) PetscCall(PetscArrayzero(e,n-1));
332: else {
333: for (i=l;i<m;i++) PetscCall(PetscArrayzero(A+l+i*ld,n-l));
334: for (i=l;i<n;i++) A[i+i*ld] = d[i];
335: }
336: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
337: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_U],&U));
338: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_V],&V));
339: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: #if !defined(PETSC_HAVE_MPIUNI)
344: static PetscErrorCode DSSynchronize_SVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
345: {
346: PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
347: PetscMPIInt n,rank,off=0,size,ldn,ld3;
348: PetscScalar *A,*U,*V;
349: PetscReal *T;
351: PetscFunctionBegin;
352: if (ds->compact) kr = 3*ld;
353: else k = (ds->n-l)*ld;
354: if (ds->state>DS_STATE_RAW) k += 2*(ds->n-l)*ld;
355: if (eigr) k += ds->n-l;
356: PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
357: PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
358: PetscCall(PetscMPIIntCast(ds->n-l,&n));
359: PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
360: PetscCall(PetscMPIIntCast(3*ld,&ld3));
361: if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
362: else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
363: if (ds->state>DS_STATE_RAW) {
364: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
365: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
366: }
367: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
368: if (!rank) {
369: if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
370: else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
371: if (ds->state>DS_STATE_RAW) {
372: PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
373: PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
374: }
375: if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
376: }
377: PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
378: if (rank) {
379: if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
380: else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
381: if (ds->state>DS_STATE_RAW) {
382: PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
383: PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
384: }
385: if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
386: }
387: if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
388: else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
389: if (ds->state>DS_STATE_RAW) {
390: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
391: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
392: }
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
395: #endif
397: static PetscErrorCode DSMatGetSize_SVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
398: {
399: DS_SVD *ctx = (DS_SVD*)ds->data;
401: PetscFunctionBegin;
402: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
403: switch (t) {
404: case DS_MAT_A:
405: *rows = ds->n;
406: *cols = ds->extrarow? ctx->m+1: ctx->m;
407: break;
408: case DS_MAT_T:
409: *rows = ds->n;
410: *cols = PetscDefined(USE_COMPLEX)? 2: 3;
411: break;
412: case DS_MAT_U:
413: *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
414: *cols = ds->n;
415: break;
416: case DS_MAT_V:
417: *rows = ds->state==DS_STATE_TRUNCATED? ctx->t: ctx->m;
418: *cols = ctx->m;
419: break;
420: default:
421: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
422: }
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: static PetscErrorCode DSSVDSetDimensions_SVD(DS ds,PetscInt m)
427: {
428: DS_SVD *ctx = (DS_SVD*)ds->data;
430: PetscFunctionBegin;
431: DSCheckAlloc(ds,1);
432: if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
433: ctx->m = ds->ld;
434: } else {
435: PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
436: ctx->m = m;
437: }
438: PetscFunctionReturn(PETSC_SUCCESS);
439: }
441: /*@
442: DSSVDSetDimensions - Sets the number of columns for a DSSVD.
444: Logically Collective
446: Input Parameters:
447: + ds - the direct solver context
448: - m - the number of columns
450: Notes:
451: This call is complementary to DSSetDimensions(), to provide a dimension
452: that is specific to this DS type.
454: Level: intermediate
456: .seealso: DSSVDGetDimensions(), DSSetDimensions()
457: @*/
458: PetscErrorCode DSSVDSetDimensions(DS ds,PetscInt m)
459: {
460: PetscFunctionBegin;
463: PetscTryMethod(ds,"DSSVDSetDimensions_C",(DS,PetscInt),(ds,m));
464: PetscFunctionReturn(PETSC_SUCCESS);
465: }
467: static PetscErrorCode DSSVDGetDimensions_SVD(DS ds,PetscInt *m)
468: {
469: DS_SVD *ctx = (DS_SVD*)ds->data;
471: PetscFunctionBegin;
472: *m = ctx->m;
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: /*@
477: DSSVDGetDimensions - Returns the number of columns for a DSSVD.
479: Not Collective
481: Input Parameter:
482: . ds - the direct solver context
484: Output Parameters:
485: . m - the number of columns
487: Level: intermediate
489: .seealso: DSSVDSetDimensions()
490: @*/
491: PetscErrorCode DSSVDGetDimensions(DS ds,PetscInt *m)
492: {
493: PetscFunctionBegin;
495: PetscAssertPointer(m,2);
496: PetscUseMethod(ds,"DSSVDGetDimensions_C",(DS,PetscInt*),(ds,m));
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: static PetscErrorCode DSDestroy_SVD(DS ds)
501: {
502: PetscFunctionBegin;
503: PetscCall(PetscFree(ds->data));
504: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSSVDSetDimensions_C",NULL));
505: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSSVDGetDimensions_C",NULL));
506: PetscFunctionReturn(PETSC_SUCCESS);
507: }
509: /*MC
510: DSSVD - Dense Singular Value Decomposition.
512: Level: beginner
514: Notes:
515: The problem is expressed as A = U*Sigma*V', where A is rectangular in
516: general, with n rows and m columns. Sigma is a diagonal matrix whose diagonal
517: elements are the arguments of DSSolve(). After solve, A is overwritten
518: with Sigma.
520: The orthogonal (or unitary) matrices of left and right singular vectors, U
521: and V, have size n and m, respectively. The number of columns m must
522: be specified via DSSVDSetDimensions().
524: If the DS object is in the intermediate state, A is assumed to be in upper
525: bidiagonal form (possibly with an arrow) and is stored in compact format
526: on matrix T. Otherwise, no particular structure is assumed. The compact
527: storage is implemented for the square case only, m=n. The extra row should
528: be interpreted in this case as an extra column.
530: Used DS matrices:
531: + DS_MAT_A - problem matrix
532: - DS_MAT_T - upper bidiagonal matrix
534: Implemented methods:
535: . 0 - Divide and Conquer (_bdsdc or _gesdd)
537: .seealso: DSCreate(), DSSetType(), DSType, DSSVDSetDimensions()
538: M*/
539: SLEPC_EXTERN PetscErrorCode DSCreate_SVD(DS ds)
540: {
541: DS_SVD *ctx;
543: PetscFunctionBegin;
544: PetscCall(PetscNew(&ctx));
545: ds->data = (void*)ctx;
547: ds->ops->allocate = DSAllocate_SVD;
548: ds->ops->view = DSView_SVD;
549: ds->ops->vectors = DSVectors_SVD;
550: ds->ops->solve[0] = DSSolve_SVD_DC;
551: ds->ops->sort = DSSort_SVD;
552: #if !defined(PETSC_HAVE_MPIUNI)
553: ds->ops->synchronize = DSSynchronize_SVD;
554: #endif
555: ds->ops->truncate = DSTruncate_SVD;
556: ds->ops->update = DSUpdateExtraRow_SVD;
557: ds->ops->destroy = DSDestroy_SVD;
558: ds->ops->matgetsize = DSMatGetSize_SVD;
559: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSSVDSetDimensions_C",DSSVDSetDimensions_SVD));
560: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSSVDGetDimensions_C",DSSVDGetDimensions_SVD));
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }