Actual source code: dssvd.c
slepc-3.18.0 2022-10-01
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: PetscErrorCode DSAllocate_SVD(DS ds,PetscInt ld)
20: {
21: DSAllocateMat_Private(ds,DS_MAT_A);
22: DSAllocateMat_Private(ds,DS_MAT_U);
23: DSAllocateMat_Private(ds,DS_MAT_V);
24: DSAllocateMat_Private(ds,DS_MAT_T);
25: PetscFree(ds->perm);
26: PetscMalloc1(ld,&ds->perm);
27: return 0;
28: }
30: /* 0 l k n-1
31: -----------------------------------------
32: |* . . |
33: | * . . |
34: | * . . |
35: | * . . |
36: | o o |
37: | o o |
38: | o o |
39: | o o |
40: | o o |
41: | o o |
42: | o x |
43: | x x |
44: | x x |
45: | x x |
46: | x x |
47: | x x |
48: | x x |
49: | x x |
50: | x x|
51: | x|
52: -----------------------------------------
53: */
55: static PetscErrorCode DSSwitchFormat_SVD(DS ds)
56: {
57: DS_SVD *ctx = (DS_SVD*)ds->data;
58: PetscReal *T;
59: PetscScalar *A;
60: PetscInt i,m=ctx->m,k=ds->k,ld=ds->ld;
63: /* switch from compact (arrow) to dense storage */
64: MatDenseGetArrayWrite(ds->omat[DS_MAT_A],&A);
65: DSGetArrayReal(ds,DS_MAT_T,&T);
66: PetscArrayzero(A,ld*ld);
67: for (i=0;i<k;i++) {
68: A[i+i*ld] = T[i];
69: A[i+k*ld] = T[i+ld];
70: }
71: A[k+k*ld] = T[k];
72: for (i=k+1;i<m;i++) {
73: A[i+i*ld] = T[i];
74: A[i-1+i*ld] = T[i-1+ld];
75: }
76: MatDenseRestoreArrayWrite(ds->omat[DS_MAT_A],&A);
77: DSRestoreArrayReal(ds,DS_MAT_T,&T);
78: return 0;
79: }
81: PetscErrorCode DSView_SVD(DS ds,PetscViewer viewer)
82: {
83: DS_SVD *ctx = (DS_SVD*)ds->data;
84: PetscViewerFormat format;
85: PetscInt i,j,r,c,m=ctx->m,rows,cols;
86: PetscReal *T,value;
88: PetscViewerGetFormat(viewer,&format);
89: if (format == PETSC_VIEWER_ASCII_INFO) return 0;
90: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
91: PetscViewerASCIIPrintf(viewer,"number of columns: %" PetscInt_FMT "\n",m);
92: return 0;
93: }
95: if (ds->compact) {
96: DSGetArrayReal(ds,DS_MAT_T,&T);
97: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
98: rows = ds->n;
99: cols = ds->extrarow? m+1: m;
100: if (format == PETSC_VIEWER_ASCII_MATLAB) {
101: PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rows,cols);
102: PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",2*ds->n);
103: PetscViewerASCIIPrintf(viewer,"zzz = [\n");
104: for (i=0;i<PetscMin(ds->n,m);i++) PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)T[i]);
105: for (i=0;i<cols-1;i++) {
106: r = PetscMax(i+2,ds->k+1);
107: c = i+1;
108: PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",c,r,(double)T[i+ds->ld]);
109: }
110: PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_T]);
111: } else {
112: for (i=0;i<rows;i++) {
113: for (j=0;j<cols;j++) {
114: if (i==j) value = T[i];
115: else if (i<ds->k && j==ds->k) value = T[PetscMin(i,j)+ds->ld];
116: else if (i+1==j && i>=ds->k) value = T[i+ds->ld];
117: else value = 0.0;
118: PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value);
119: }
120: PetscViewerASCIIPrintf(viewer,"\n");
121: }
122: }
123: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
124: PetscViewerFlush(viewer);
125: DSRestoreArrayReal(ds,DS_MAT_T,&T);
126: } else DSViewMat(ds,viewer,DS_MAT_A);
127: if (ds->state>DS_STATE_INTERMEDIATE) {
128: DSViewMat(ds,viewer,DS_MAT_U);
129: DSViewMat(ds,viewer,DS_MAT_V);
130: }
131: return 0;
132: }
134: PetscErrorCode DSVectors_SVD(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
135: {
136: switch (mat) {
137: case DS_MAT_U:
138: case DS_MAT_V:
139: if (rnorm) *rnorm = 0.0;
140: break;
141: default:
142: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
143: }
144: return 0;
145: }
147: PetscErrorCode DSSort_SVD(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
148: {
149: DS_SVD *ctx = (DS_SVD*)ds->data;
150: PetscInt n,l,i,*perm,ld=ds->ld;
151: PetscScalar *A;
152: PetscReal *d;
154: if (!ds->sc) return 0;
156: l = ds->l;
157: n = PetscMin(ds->n,ctx->m);
158: DSGetArrayReal(ds,DS_MAT_T,&d);
159: perm = ds->perm;
160: if (!rr) DSSortEigenvaluesReal_Private(ds,d,perm);
161: else DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE);
162: for (i=l;i<n;i++) wr[i] = d[perm[i]];
163: DSPermuteBoth_Private(ds,l,n,ds->n,ctx->m,DS_MAT_U,DS_MAT_V,perm);
164: for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
165: if (!ds->compact) {
166: MatDenseGetArray(ds->omat[DS_MAT_A],&A);
167: for (i=l;i<n;i++) A[i+i*ld] = wr[i];
168: MatDenseRestoreArray(ds->omat[DS_MAT_A],&A);
169: }
170: DSRestoreArrayReal(ds,DS_MAT_T,&d);
171: return 0;
172: }
174: PetscErrorCode DSUpdateExtraRow_SVD(DS ds)
175: {
176: DS_SVD *ctx = (DS_SVD*)ds->data;
177: PetscInt i;
178: PetscBLASInt n=0,m=0,ld,incx=1;
179: PetscScalar *A,*x,*y,one=1.0,zero=0.0;
180: PetscReal *T,*e,beta;
181: const PetscScalar *U;
184: PetscBLASIntCast(ds->n,&n);
185: PetscBLASIntCast(ctx->m,&m);
186: PetscBLASIntCast(ds->ld,&ld);
187: MatDenseGetArrayRead(ds->omat[DS_MAT_U],&U);
188: if (ds->compact) {
189: DSGetArrayReal(ds,DS_MAT_T,&T);
190: e = T+ld;
191: beta = e[m-1]; /* in compact, we assume all entries are zero except the last one */
192: for (i=0;i<n;i++) e[i] = PetscRealPart(beta*U[n-1+i*ld]);
193: ds->k = m;
194: DSRestoreArrayReal(ds,DS_MAT_T,&T);
195: } else {
196: MatDenseGetArray(ds->omat[DS_MAT_A],&A);
197: DSAllocateWork_Private(ds,2*ld,0,0);
198: x = ds->work;
199: y = ds->work+ld;
200: for (i=0;i<n;i++) x[i] = PetscConj(A[i+m*ld]);
201: PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,U,&ld,x,&incx,&zero,y,&incx));
202: for (i=0;i<n;i++) A[i+m*ld] = PetscConj(y[i]);
203: ds->k = m;
204: MatDenseRestoreArray(ds->omat[DS_MAT_A],&A);
205: }
206: MatDenseRestoreArrayRead(ds->omat[DS_MAT_U],&U);
207: return 0;
208: }
210: PetscErrorCode DSTruncate_SVD(DS ds,PetscInt n,PetscBool trim)
211: {
212: PetscInt i,ld=ds->ld,l=ds->l;
213: PetscScalar *A;
214: DS_SVD *ctx = (DS_SVD*)ds->data;
216: if (!ds->compact && ds->extrarow) MatDenseGetArray(ds->omat[DS_MAT_A],&A);
217: if (trim) {
218: if (!ds->compact && ds->extrarow) { /* clean extra column */
219: for (i=l;i<ds->n;i++) A[i+ctx->m*ld] = 0.0;
220: }
221: ds->l = 0;
222: ds->k = 0;
223: ds->n = n;
224: ctx->m = n;
225: ds->t = ds->n; /* truncated length equal to the new dimension */
226: ctx->t = ctx->m; /* must also keep the previous dimension of V */
227: } else {
228: if (!ds->compact && ds->extrarow && ds->k==ds->n) {
229: /* copy entries of extra column to the new position, then clean last row */
230: for (i=l;i<n;i++) A[i+n*ld] = A[i+ctx->m*ld];
231: for (i=l;i<ds->n;i++) A[i+ctx->m*ld] = 0.0;
232: }
233: ds->k = (ds->extrarow)? n: 0;
234: ds->t = ds->n; /* truncated length equal to previous dimension */
235: ctx->t = ctx->m; /* must also keep the previous dimension of V */
236: ds->n = n;
237: ctx->m = n;
238: }
239: if (!ds->compact && ds->extrarow) MatDenseRestoreArray(ds->omat[DS_MAT_A],&A);
240: return 0;
241: }
243: PetscErrorCode DSSolve_SVD_DC(DS ds,PetscScalar *wr,PetscScalar *wi)
244: {
245: DS_SVD *ctx = (DS_SVD*)ds->data;
246: PetscInt i,j;
247: PetscBLASInt n1,m1,info,l = 0,n = 0,m = 0,nm,ld,off,lwork;
248: PetscScalar *A,*U,*V,*W,qwork;
249: PetscReal *d,*e,*Ur,*Vr;
252: PetscBLASIntCast(ds->n,&n);
253: PetscBLASIntCast(ctx->m,&m);
254: PetscBLASIntCast(ds->l,&l);
255: PetscBLASIntCast(ds->ld,&ld);
256: n1 = n-l; /* n1 = size of leading block, excl. locked + size of trailing block */
257: m1 = m-l;
258: off = l+l*ld;
259: MatDenseGetArray(ds->omat[DS_MAT_A],&A);
260: MatDenseGetArrayWrite(ds->omat[DS_MAT_U],&U);
261: MatDenseGetArrayWrite(ds->omat[DS_MAT_V],&V);
262: DSGetArrayReal(ds,DS_MAT_T,&d);
263: e = d+ld;
264: PetscArrayzero(U,ld*ld);
265: for (i=0;i<l;i++) U[i+i*ld] = 1.0;
266: PetscArrayzero(V,ld*ld);
267: for (i=0;i<l;i++) V[i+i*ld] = 1.0;
269: if (ds->state>DS_STATE_RAW) {
270: /* solve bidiagonal SVD problem */
271: for (i=0;i<l;i++) wr[i] = d[i];
272: #if defined(PETSC_USE_COMPLEX)
273: DSAllocateWork_Private(ds,0,3*n1*n1+4*n1+2*ld*ld,8*n1);
274: Ur = ds->rwork+3*n1*n1+4*n1;
275: Vr = ds->rwork+3*n1*n1+4*n1+ld*ld;
276: #else
277: DSAllocateWork_Private(ds,0,3*n1*n1+4*n1+ld*ld,8*n1);
278: Ur = U;
279: Vr = ds->rwork+3*n1*n1+4*n1;
280: #endif
281: PetscCallBLAS("LAPACKbdsdc",LAPACKbdsdc_("U","I",&n1,d+l,e+l,Ur+off,&ld,Vr+off,&ld,NULL,NULL,ds->rwork,ds->iwork,&info));
282: SlepcCheckLapackInfo("bdsdc",info);
283: for (i=l;i<n;i++) {
284: for (j=l;j<n;j++) {
285: #if defined(PETSC_USE_COMPLEX)
286: U[i+j*ld] = Ur[i+j*ld];
287: #endif
288: V[i+j*ld] = PetscConj(Vr[j+i*ld]); /* transpose VT returned by Lapack */
289: }
290: }
291: } else {
292: /* solve general rectangular SVD problem */
293: DSAllocateMat_Private(ds,DS_MAT_W);
294: MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&W);
295: if (ds->compact) DSSwitchFormat_SVD(ds);
296: for (i=0;i<l;i++) wr[i] = d[i];
297: nm = PetscMin(n,m);
298: DSAllocateWork_Private(ds,0,0,8*nm);
299: lwork = -1;
300: #if defined(PETSC_USE_COMPLEX)
301: DSAllocateWork_Private(ds,0,5*nm*nm+7*nm,0);
302: PetscCallBLAS("LAPACKgesdd",LAPACKgesdd_("A",&n1,&m1,A+off,&ld,d+l,U+off,&ld,W+off,&ld,&qwork,&lwork,ds->rwork,ds->iwork,&info));
303: #else
304: PetscCallBLAS("LAPACKgesdd",LAPACKgesdd_("A",&n1,&m1,A+off,&ld,d+l,U+off,&ld,W+off,&ld,&qwork,&lwork,ds->iwork,&info));
305: #endif
306: SlepcCheckLapackInfo("gesdd",info);
307: PetscBLASIntCast((PetscInt)PetscRealPart(qwork),&lwork);
308: DSAllocateWork_Private(ds,lwork,0,0);
309: #if defined(PETSC_USE_COMPLEX)
310: 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));
311: #else
312: PetscCallBLAS("LAPACKgesdd",LAPACKgesdd_("A",&n1,&m1,A+off,&ld,d+l,U+off,&ld,W+off,&ld,ds->work,&lwork,ds->iwork,&info));
313: #endif
314: SlepcCheckLapackInfo("gesdd",info);
315: for (i=l;i<m;i++) {
316: for (j=l;j<m;j++) V[i+j*ld] = PetscConj(W[j+i*ld]); /* transpose VT returned by Lapack */
317: }
318: MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&W);
319: }
320: for (i=l;i<PetscMin(ds->n,ctx->m);i++) wr[i] = d[i];
322: /* create diagonal matrix as a result */
323: if (ds->compact) PetscArrayzero(e,n-1);
324: else {
325: for (i=l;i<m;i++) PetscArrayzero(A+l+i*ld,n-l);
326: for (i=l;i<n;i++) A[i+i*ld] = d[i];
327: }
328: MatDenseRestoreArray(ds->omat[DS_MAT_A],&A);
329: MatDenseRestoreArrayWrite(ds->omat[DS_MAT_U],&U);
330: MatDenseRestoreArrayWrite(ds->omat[DS_MAT_V],&V);
331: DSRestoreArrayReal(ds,DS_MAT_T,&d);
332: return 0;
333: }
335: PetscErrorCode DSSynchronize_SVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
336: {
337: PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
338: PetscMPIInt n,rank,off=0,size,ldn,ld3;
339: PetscScalar *A,*U,*V;
340: PetscReal *T;
342: if (ds->compact) kr = 3*ld;
343: else k = (ds->n-l)*ld;
344: if (ds->state>DS_STATE_RAW) k += 2*(ds->n-l)*ld;
345: if (eigr) k += ds->n-l;
346: DSAllocateWork_Private(ds,k+kr,0,0);
347: PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size);
348: PetscMPIIntCast(ds->n-l,&n);
349: PetscMPIIntCast(ld*(ds->n-l),&ldn);
350: PetscMPIIntCast(3*ld,&ld3);
351: if (ds->compact) DSGetArrayReal(ds,DS_MAT_T,&T);
352: else MatDenseGetArray(ds->omat[DS_MAT_A],&A);
353: if (ds->state>DS_STATE_RAW) {
354: MatDenseGetArray(ds->omat[DS_MAT_U],&U);
355: MatDenseGetArray(ds->omat[DS_MAT_V],&V);
356: }
357: MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
358: if (!rank) {
359: if (ds->compact) MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
360: else MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
361: if (ds->state>DS_STATE_RAW) {
362: MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
363: MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
364: }
365: if (eigr) MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
366: }
367: MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
368: if (rank) {
369: if (ds->compact) MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds));
370: else MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
371: if (ds->state>DS_STATE_RAW) {
372: MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
373: MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
374: }
375: if (eigr) MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
376: }
377: if (ds->compact) DSRestoreArrayReal(ds,DS_MAT_T,&T);
378: MatDenseRestoreArray(ds->omat[DS_MAT_A],&A);
379: if (ds->state>DS_STATE_RAW) {
380: MatDenseRestoreArray(ds->omat[DS_MAT_U],&U);
381: MatDenseRestoreArray(ds->omat[DS_MAT_V],&V);
382: }
383: return 0;
384: }
386: PetscErrorCode DSMatGetSize_SVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
387: {
388: DS_SVD *ctx = (DS_SVD*)ds->data;
391: switch (t) {
392: case DS_MAT_A:
393: *rows = ds->n;
394: *cols = ds->extrarow? ctx->m+1: ctx->m;
395: break;
396: case DS_MAT_T:
397: *rows = ds->n;
398: *cols = PetscDefined(USE_COMPLEX)? 2: 3;
399: break;
400: case DS_MAT_U:
401: *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
402: *cols = ds->n;
403: break;
404: case DS_MAT_V:
405: *rows = ds->state==DS_STATE_TRUNCATED? ctx->t: ctx->m;
406: *cols = ctx->m;
407: break;
408: default:
409: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
410: }
411: return 0;
412: }
414: static PetscErrorCode DSSVDSetDimensions_SVD(DS ds,PetscInt m)
415: {
416: DS_SVD *ctx = (DS_SVD*)ds->data;
418: DSCheckAlloc(ds,1);
419: if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
420: ctx->m = ds->ld;
421: } else {
423: ctx->m = m;
424: }
425: return 0;
426: }
428: /*@
429: DSSVDSetDimensions - Sets the number of columns for a DSSVD.
431: Logically Collective on ds
433: Input Parameters:
434: + ds - the direct solver context
435: - m - the number of columns
437: Notes:
438: This call is complementary to DSSetDimensions(), to provide a dimension
439: that is specific to this DS type.
441: Level: intermediate
443: .seealso: DSSVDGetDimensions(), DSSetDimensions()
444: @*/
445: PetscErrorCode DSSVDSetDimensions(DS ds,PetscInt m)
446: {
449: PetscTryMethod(ds,"DSSVDSetDimensions_C",(DS,PetscInt),(ds,m));
450: return 0;
451: }
453: static PetscErrorCode DSSVDGetDimensions_SVD(DS ds,PetscInt *m)
454: {
455: DS_SVD *ctx = (DS_SVD*)ds->data;
457: *m = ctx->m;
458: return 0;
459: }
461: /*@
462: DSSVDGetDimensions - Returns the number of columns for a DSSVD.
464: Not collective
466: Input Parameter:
467: . ds - the direct solver context
469: Output Parameters:
470: . m - the number of columns
472: Level: intermediate
474: .seealso: DSSVDSetDimensions()
475: @*/
476: PetscErrorCode DSSVDGetDimensions(DS ds,PetscInt *m)
477: {
480: PetscUseMethod(ds,"DSSVDGetDimensions_C",(DS,PetscInt*),(ds,m));
481: return 0;
482: }
484: PetscErrorCode DSDestroy_SVD(DS ds)
485: {
486: PetscFree(ds->data);
487: PetscObjectComposeFunction((PetscObject)ds,"DSSVDSetDimensions_C",NULL);
488: PetscObjectComposeFunction((PetscObject)ds,"DSSVDGetDimensions_C",NULL);
489: return 0;
490: }
492: /*MC
493: DSSVD - Dense Singular Value Decomposition.
495: Level: beginner
497: Notes:
498: The problem is expressed as A = U*Sigma*V', where A is rectangular in
499: general, with n rows and m columns. Sigma is a diagonal matrix whose diagonal
500: elements are the arguments of DSSolve(). After solve, A is overwritten
501: with Sigma.
503: The orthogonal (or unitary) matrices of left and right singular vectors, U
504: and V, have size n and m, respectively. The number of columns m must
505: be specified via DSSVDSetDimensions().
507: If the DS object is in the intermediate state, A is assumed to be in upper
508: bidiagonal form (possibly with an arrow) and is stored in compact format
509: on matrix T. Otherwise, no particular structure is assumed. The compact
510: storage is implemented for the square case only, m=n. The extra row should
511: be interpreted in this case as an extra column.
513: Used DS matrices:
514: + DS_MAT_A - problem matrix
515: - DS_MAT_T - upper bidiagonal matrix
517: Implemented methods:
518: . 0 - Divide and Conquer (_bdsdc or _gesdd)
520: .seealso: DSCreate(), DSSetType(), DSType, DSSVDSetDimensions()
521: M*/
522: SLEPC_EXTERN PetscErrorCode DSCreate_SVD(DS ds)
523: {
524: DS_SVD *ctx;
526: PetscNew(&ctx);
527: ds->data = (void*)ctx;
529: ds->ops->allocate = DSAllocate_SVD;
530: ds->ops->view = DSView_SVD;
531: ds->ops->vectors = DSVectors_SVD;
532: ds->ops->solve[0] = DSSolve_SVD_DC;
533: ds->ops->sort = DSSort_SVD;
534: ds->ops->synchronize = DSSynchronize_SVD;
535: ds->ops->truncate = DSTruncate_SVD;
536: ds->ops->update = DSUpdateExtraRow_SVD;
537: ds->ops->destroy = DSDestroy_SVD;
538: ds->ops->matgetsize = DSMatGetSize_SVD;
539: PetscObjectComposeFunction((PetscObject)ds,"DSSVDSetDimensions_C",DSSVDSetDimensions_SVD);
540: PetscObjectComposeFunction((PetscObject)ds,"DSSVDGetDimensions_C",DSSVDGetDimensions_SVD);
541: return 0;
542: }