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             : } DS_SVD;
      18             : 
      19          76 : static PetscErrorCode DSAllocate_SVD(DS ds,PetscInt ld)
      20             : {
      21          76 :   PetscFunctionBegin;
      22          76 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
      23          76 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
      24          76 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_V));
      25          76 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
      26          76 :   PetscCall(PetscFree(ds->perm));
      27          76 :   PetscCall(PetscMalloc1(ld,&ds->perm));
      28          76 :   PetscFunctionReturn(PETSC_SUCCESS);
      29             : }
      30             : 
      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             : */
      55             : 
      56         471 : static PetscErrorCode DSSwitchFormat_SVD(DS ds)
      57             : {
      58         471 :   DS_SVD         *ctx = (DS_SVD*)ds->data;
      59         471 :   PetscReal      *T;
      60         471 :   PetscScalar    *A;
      61         471 :   PetscInt       i,m=ctx->m,k=ds->k,ld=ds->ld;
      62             : 
      63         471 :   PetscFunctionBegin;
      64         471 :   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         471 :   PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_A],&A));
      67         471 :   PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
      68         471 :   PetscCall(PetscArrayzero(A,ld*ld));
      69        2714 :   for (i=0;i<k;i++) {
      70        2243 :     A[i+i*ld] = T[i];
      71        2243 :     A[i+k*ld] = T[i+ld];
      72             :   }
      73         471 :   A[k+k*ld] = T[k];
      74        2373 :   for (i=k+1;i<m;i++) {
      75        1902 :     A[i+i*ld]   = T[i];
      76        1902 :     A[i-1+i*ld] = T[i-1+ld];
      77             :   }
      78         471 :   PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_A],&A));
      79         471 :   PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
      80         471 :   PetscFunctionReturn(PETSC_SUCCESS);
      81             : }
      82             : 
      83          11 : static PetscErrorCode DSView_SVD(DS ds,PetscViewer viewer)
      84             : {
      85          11 :   DS_SVD            *ctx = (DS_SVD*)ds->data;
      86          11 :   PetscViewerFormat format;
      87          11 :   PetscInt          i,j,r,c,m=ctx->m,rows,cols;
      88          11 :   PetscReal         *T,value;
      89             : 
      90          11 :   PetscFunctionBegin;
      91          11 :   PetscCall(PetscViewerGetFormat(viewer,&format));
      92          11 :   if (format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
      93          10 :   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
      94           5 :     PetscCall(PetscViewerASCIIPrintf(viewer,"number of columns: %" PetscInt_FMT "\n",m));
      95           5 :     PetscFunctionReturn(PETSC_SUCCESS);
      96             :   }
      97           5 :   PetscCheck(m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
      98           5 :   if (ds->compact) {
      99           3 :     PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
     100           3 :     PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
     101           3 :     rows = ds->n;
     102           3 :     cols = ds->extrarow? m+1: m;
     103           3 :     if (format == PETSC_VIEWER_ASCII_MATLAB) {
     104           3 :       PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rows,cols));
     105           3 :       PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",2*ds->n));
     106           3 :       PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
     107          33 :       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          31 :       for (i=0;i<cols-1;i++) {
     109          28 :         r = PetscMax(i+2,ds->k+1);
     110          28 :         c = i+1;
     111          28 :         PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT "  %18.16e\n",c,r,(double)T[i+ds->ld]));
     112             :       }
     113           3 :       PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_T]));
     114             :     } else {
     115           0 :       for (i=0;i<rows;i++) {
     116           0 :         for (j=0;j<cols;j++) {
     117           0 :           if (i==j) value = T[i];
     118           0 :           else if (i<ds->k && j==ds->k) value = T[PetscMin(i,j)+ds->ld];
     119           0 :           else if (i+1==j && i>=ds->k) value = T[i+ds->ld];
     120             :           else value = 0.0;
     121           0 :           PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
     122             :         }
     123           0 :         PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
     124             :       }
     125             :     }
     126           3 :     PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
     127           3 :     PetscCall(PetscViewerFlush(viewer));
     128           3 :     PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
     129           2 :   } else PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
     130           5 :   if (ds->state>DS_STATE_INTERMEDIATE) {
     131           0 :     PetscCall(DSViewMat(ds,viewer,DS_MAT_U));
     132           0 :     PetscCall(DSViewMat(ds,viewer,DS_MAT_V));
     133             :   }
     134           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     135             : }
     136             : 
     137           2 : static PetscErrorCode DSVectors_SVD(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
     138             : {
     139           2 :   PetscFunctionBegin;
     140           2 :   switch (mat) {
     141           2 :     case DS_MAT_U:
     142             :     case DS_MAT_V:
     143           2 :       if (rnorm) *rnorm = 0.0;
     144           2 :       break;
     145           0 :     default:
     146           0 :       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
     147             :   }
     148           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     149             : }
     150             : 
     151        1372 : static PetscErrorCode DSSort_SVD(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
     152             : {
     153        1372 :   DS_SVD         *ctx = (DS_SVD*)ds->data;
     154        1372 :   PetscInt       n,l,i,*perm,ld=ds->ld;
     155        1372 :   PetscScalar    *A;
     156        1372 :   PetscReal      *d;
     157             : 
     158        1372 :   PetscFunctionBegin;
     159        1372 :   if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
     160        1372 :   PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
     161        1372 :   l = ds->l;
     162        1372 :   n = PetscMin(ds->n,ctx->m);
     163        1372 :   PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
     164        1372 :   perm = ds->perm;
     165        1372 :   if (!rr) PetscCall(DSSortEigenvaluesReal_Private(ds,d,perm));
     166           0 :   else PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
     167       14863 :   for (i=l;i<n;i++) wr[i] = d[perm[i]];
     168        1372 :   PetscCall(DSPermuteBoth_Private(ds,l,n,ds->n,ctx->m,DS_MAT_U,DS_MAT_V,perm));
     169       14863 :   for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
     170        1372 :   if (!ds->compact) {
     171         290 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     172        3526 :     for (i=l;i<n;i++) A[i+i*ld] = wr[i];
     173         290 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     174             :   }
     175        1372 :   PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
     176        1372 :   PetscFunctionReturn(PETSC_SUCCESS);
     177             : }
     178             : 
     179        1081 : static PetscErrorCode DSUpdateExtraRow_SVD(DS ds)
     180             : {
     181        1081 :   DS_SVD            *ctx = (DS_SVD*)ds->data;
     182        1081 :   PetscInt          i;
     183        1081 :   PetscBLASInt      n=0,m=0,ld,incx=1;
     184        1081 :   PetscScalar       *A,*x,*y,one=1.0,zero=0.0;
     185        1081 :   PetscReal         *T,*e,beta;
     186        1081 :   const PetscScalar *U;
     187             : 
     188        1081 :   PetscFunctionBegin;
     189        1081 :   PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
     190        1081 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     191        1081 :   PetscCall(PetscBLASIntCast(ctx->m,&m));
     192        1081 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     193        1081 :   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_U],&U));
     194        1081 :   if (ds->compact) {
     195        1080 :     PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
     196        1080 :     e = T+ld;
     197        1080 :     beta = e[m-1];   /* in compact, we assume all entries are zero except the last one */
     198       11868 :     for (i=0;i<n;i++) e[i] = PetscRealPart(beta*U[n-1+i*ld]);
     199        1080 :     ds->k = m;
     200        1080 :     PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
     201             :   } else {
     202           1 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     203           1 :     PetscCall(DSAllocateWork_Private(ds,2*ld,0,0));
     204           1 :     x = ds->work;
     205           1 :     y = ds->work+ld;
     206          16 :     for (i=0;i<n;i++) x[i] = PetscConj(A[i+m*ld]);
     207           1 :     PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,U,&ld,x,&incx,&zero,y,&incx));
     208          16 :     for (i=0;i<n;i++) A[i+m*ld] = PetscConj(y[i]);
     209           1 :     ds->k = m;
     210           1 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     211             :   }
     212        1081 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_U],&U));
     213        1081 :   PetscFunctionReturn(PETSC_SUCCESS);
     214             : }
     215             : 
     216         508 : static PetscErrorCode DSTruncate_SVD(DS ds,PetscInt n,PetscBool trim)
     217             : {
     218         508 :   PetscInt    i,ld=ds->ld,l=ds->l;
     219         508 :   PetscScalar *A;
     220         508 :   DS_SVD      *ctx = (DS_SVD*)ds->data;
     221             : 
     222         508 :   PetscFunctionBegin;
     223         508 :   if (!ds->compact && ds->extrarow) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     224         508 :   if (trim) {
     225          39 :     if (!ds->compact && ds->extrarow) {   /* clean extra column */
     226           0 :       for (i=l;i<ds->n;i++) A[i+ctx->m*ld] = 0.0;
     227             :     }
     228          39 :     ds->l  = 0;
     229          39 :     ds->k  = 0;
     230          39 :     ds->n  = n;
     231          39 :     ctx->m = n;
     232          39 :     ds->t  = ds->n;   /* truncated length equal to the new dimension */
     233          39 :     ctx->t = ctx->m;  /* must also keep the previous dimension of V */
     234             :   } else {
     235         469 :     if (!ds->compact && ds->extrarow && ds->k==ds->n) {
     236             :       /* copy entries of extra column to the new position, then clean last row */
     237           0 :       for (i=l;i<n;i++) A[i+n*ld] = A[i+ctx->m*ld];
     238           0 :       for (i=l;i<ds->n;i++) A[i+ctx->m*ld] = 0.0;
     239             :     }
     240         469 :     ds->k  = (ds->extrarow)? n: 0;
     241         469 :     ds->t  = ds->n;   /* truncated length equal to previous dimension */
     242         469 :     ctx->t = ctx->m;  /* must also keep the previous dimension of V */
     243         469 :     ds->n  = n;
     244         469 :     ctx->m = n;
     245             :   }
     246         508 :   if (!ds->compact && ds->extrarow) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     247         508 :   PetscFunctionReturn(PETSC_SUCCESS);
     248             : }
     249             : 
     250        1428 : static PetscErrorCode DSSolve_SVD_DC(DS ds,PetscScalar *wr,PetscScalar *wi)
     251             : {
     252        1428 :   DS_SVD         *ctx = (DS_SVD*)ds->data;
     253        1428 :   PetscInt       i,j;
     254        1428 :   PetscBLASInt   n1,m1,info,l = 0,n = 0,m = 0,nm,ld,off,lwork;
     255        1428 :   PetscScalar    *A,*U,*V,*W,qwork;
     256        1428 :   PetscReal      *d,*e,*Ur,*Vr;
     257             : 
     258        1428 :   PetscFunctionBegin;
     259        1428 :   PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
     260        1428 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     261        1428 :   PetscCall(PetscBLASIntCast(ctx->m,&m));
     262        1428 :   PetscCall(PetscBLASIntCast(ds->l,&l));
     263        1428 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     264        1428 :   n1 = n-l;     /* n1 = size of leading block, excl. locked + size of trailing block */
     265        1428 :   m1 = m-l;
     266        1428 :   off = l+l*ld;
     267        1428 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     268        1428 :   PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_U],&U));
     269        1428 :   PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_V],&V));
     270        1428 :   PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
     271        1428 :   e = d+ld;
     272        1428 :   PetscCall(PetscArrayzero(U,ld*ld));
     273        2013 :   for (i=0;i<l;i++) U[i+i*ld] = 1.0;
     274        1428 :   PetscCall(PetscArrayzero(V,ld*ld));
     275        2013 :   for (i=0;i<l;i++) V[i+i*ld] = 1.0;
     276             : 
     277        1428 :   if (ds->state>DS_STATE_RAW) {
     278             :     /* solve bidiagonal SVD problem */
     279        1046 :     for (i=0;i<l;i++) wr[i] = d[i];
     280             : #if defined(PETSC_USE_COMPLEX)
     281         611 :     PetscCall(DSAllocateWork_Private(ds,0,3*n1*n1+4*n1+2*ld*ld,8*n1));
     282         611 :     Ur = ds->rwork+3*n1*n1+4*n1;
     283         611 :     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         611 :     PetscCallBLAS("LAPACKbdsdc",LAPACKbdsdc_("U","I",&n1,d+l,e+l,Ur+off,&ld,Vr+off,&ld,NULL,NULL,ds->rwork,ds->iwork,&info));
     290         611 :     SlepcCheckLapackInfo("bdsdc",info);
     291        6368 :     for (i=l;i<n;i++) {
     292       60894 :       for (j=l;j<n;j++) {
     293             : #if defined(PETSC_USE_COMPLEX)
     294       55137 :         U[i+j*ld] = Ur[i+j*ld];
     295             : #endif
     296       55137 :         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         817 :     PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
     302         817 :     PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&W));
     303         817 :     if (ds->compact) PetscCall(DSSwitchFormat_SVD(ds));
     304         967 :     for (i=0;i<l;i++) wr[i] = d[i];
     305         817 :     nm = PetscMin(n,m);
     306         817 :     PetscCall(DSAllocateWork_Private(ds,0,0,8*nm));
     307         817 :     lwork = -1;
     308             : #if defined(PETSC_USE_COMPLEX)
     309         817 :     PetscCall(DSAllocateWork_Private(ds,0,5*nm*nm+7*nm,0));
     310         817 :     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         817 :     SlepcCheckLapackInfo("gesdd",info);
     315         817 :     PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(qwork),&lwork));
     316         817 :     PetscCall(DSAllocateWork_Private(ds,lwork,0,0));
     317             : #if defined(PETSC_USE_COMPLEX)
     318         817 :     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         817 :     SlepcCheckLapackInfo("gesdd",info);
     323        9521 :     for (i=l;i<m;i++) {
     324      116058 :       for (j=l;j<m;j++) V[i+j*ld] = PetscConj(W[j+i*ld]);  /* transpose VT returned by Lapack */
     325             :     }
     326         817 :     PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&W));
     327             :   }
     328       15883 :   for (i=l;i<PetscMin(ds->n,ctx->m);i++) wr[i] = d[i];
     329             : 
     330             :   /* create diagonal matrix as a result */
     331        1428 :   if (ds->compact) PetscCall(PetscArrayzero(e,n-1));
     332             :   else {
     333        4552 :     for (i=l;i<m;i++) PetscCall(PetscArrayzero(A+l+i*ld,n-l));
     334        4580 :     for (i=l;i<n;i++) A[i+i*ld] = d[i];
     335             :   }
     336        1428 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     337        1428 :   PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_U],&U));
     338        1428 :   PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_V],&V));
     339        1428 :   PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
     340        1428 :   PetscFunctionReturn(PETSC_SUCCESS);
     341             : }
     342             : 
     343             : #if !defined(PETSC_HAVE_MPIUNI)
     344           6 : static PetscErrorCode DSSynchronize_SVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
     345             : {
     346           6 :   PetscInt       ld=ds->ld,l=ds->l,k=0,kr=0;
     347           6 :   PetscMPIInt    n,rank,off=0,size,ldn,ld3;
     348           6 :   PetscScalar    *A,*U,*V;
     349           6 :   PetscReal      *T;
     350             : 
     351           6 :   PetscFunctionBegin;
     352           6 :   if (ds->compact) kr = 3*ld;
     353           0 :   else k = (ds->n-l)*ld;
     354           6 :   if (ds->state>DS_STATE_RAW) k += 2*(ds->n-l)*ld;
     355           6 :   if (eigr) k += ds->n-l;
     356           6 :   PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
     357           6 :   PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
     358           6 :   PetscCall(PetscMPIIntCast(ds->n-l,&n));
     359           6 :   PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
     360           6 :   PetscCall(PetscMPIIntCast(3*ld,&ld3));
     361           6 :   if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
     362           0 :   else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     363           6 :   if (ds->state>DS_STATE_RAW) {
     364           6 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
     365           6 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
     366             :   }
     367           6 :   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
     368           6 :   if (!rank) {
     369           3 :     if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     370           0 :     else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     371           3 :     if (ds->state>DS_STATE_RAW) {
     372           3 :       PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     373           3 :       PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     374             :     }
     375           3 :     if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     376             :   }
     377          12 :   PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
     378           6 :   if (rank) {
     379           3 :     if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
     380           0 :     else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     381           3 :     if (ds->state>DS_STATE_RAW) {
     382           3 :       PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     383           3 :       PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     384             :     }
     385           3 :     if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     386             :   }
     387           6 :   if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
     388           0 :   else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     389           6 :   if (ds->state>DS_STATE_RAW) {
     390           6 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
     391           6 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
     392             :   }
     393           6 :   PetscFunctionReturn(PETSC_SUCCESS);
     394             : }
     395             : #endif
     396             : 
     397        3050 : static PetscErrorCode DSMatGetSize_SVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
     398             : {
     399        3050 :   DS_SVD *ctx = (DS_SVD*)ds->data;
     400             : 
     401        3050 :   PetscFunctionBegin;
     402        3050 :   PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSSVDSetDimensions()");
     403        3050 :   switch (t) {
     404         346 :     case DS_MAT_A:
     405         346 :       *rows = ds->n;
     406         346 :       *cols = ds->extrarow? ctx->m+1: ctx->m;
     407         346 :       break;
     408           0 :     case DS_MAT_T:
     409           0 :       *rows = ds->n;
     410           0 :       *cols = PetscDefined(USE_COMPLEX)? 2: 3;
     411           0 :       break;
     412        1353 :     case DS_MAT_U:
     413        1353 :       *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
     414        1353 :       *cols = ds->n;
     415        1353 :       break;
     416        1351 :     case DS_MAT_V:
     417        1351 :       *rows = ds->state==DS_STATE_TRUNCATED? ctx->t: ctx->m;
     418        1351 :       *cols = ctx->m;
     419        1351 :       break;
     420           0 :     default:
     421           0 :       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
     422             :   }
     423        3050 :   PetscFunctionReturn(PETSC_SUCCESS);
     424             : }
     425             : 
     426        1428 : static PetscErrorCode DSSVDSetDimensions_SVD(DS ds,PetscInt m)
     427             : {
     428        1428 :   DS_SVD *ctx = (DS_SVD*)ds->data;
     429             : 
     430        1428 :   PetscFunctionBegin;
     431        1428 :   DSCheckAlloc(ds,1);
     432        1428 :   if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
     433           0 :     ctx->m = ds->ld;
     434             :   } else {
     435        1428 :     PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
     436        1428 :     ctx->m = m;
     437             :   }
     438        1428 :   PetscFunctionReturn(PETSC_SUCCESS);
     439             : }
     440             : 
     441             : /*@
     442             :    DSSVDSetDimensions - Sets the number of columns for a DSSVD.
     443             : 
     444             :    Logically Collective
     445             : 
     446             :    Input Parameters:
     447             : +  ds - the direct solver context
     448             : -  m  - the number of columns
     449             : 
     450             :    Notes:
     451             :    This call is complementary to DSSetDimensions(), to provide a dimension
     452             :    that is specific to this DS type.
     453             : 
     454             :    Level: intermediate
     455             : 
     456             : .seealso: DSSVDGetDimensions(), DSSetDimensions()
     457             : @*/
     458        1428 : PetscErrorCode DSSVDSetDimensions(DS ds,PetscInt m)
     459             : {
     460        1428 :   PetscFunctionBegin;
     461        1428 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     462        5712 :   PetscValidLogicalCollectiveInt(ds,m,2);
     463        1428 :   PetscTryMethod(ds,"DSSVDSetDimensions_C",(DS,PetscInt),(ds,m));
     464        1428 :   PetscFunctionReturn(PETSC_SUCCESS);
     465             : }
     466             : 
     467           2 : static PetscErrorCode DSSVDGetDimensions_SVD(DS ds,PetscInt *m)
     468             : {
     469           2 :   DS_SVD *ctx = (DS_SVD*)ds->data;
     470             : 
     471           2 :   PetscFunctionBegin;
     472           2 :   *m = ctx->m;
     473           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     474             : }
     475             : 
     476             : /*@
     477             :    DSSVDGetDimensions - Returns the number of columns for a DSSVD.
     478             : 
     479             :    Not Collective
     480             : 
     481             :    Input Parameter:
     482             : .  ds - the direct solver context
     483             : 
     484             :    Output Parameters:
     485             : .  m - the number of columns
     486             : 
     487             :    Level: intermediate
     488             : 
     489             : .seealso: DSSVDSetDimensions()
     490             : @*/
     491           2 : PetscErrorCode DSSVDGetDimensions(DS ds,PetscInt *m)
     492             : {
     493           2 :   PetscFunctionBegin;
     494           2 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     495           2 :   PetscAssertPointer(m,2);
     496           2 :   PetscUseMethod(ds,"DSSVDGetDimensions_C",(DS,PetscInt*),(ds,m));
     497           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     498             : }
     499             : 
     500          73 : static PetscErrorCode DSDestroy_SVD(DS ds)
     501             : {
     502          73 :   PetscFunctionBegin;
     503          73 :   PetscCall(PetscFree(ds->data));
     504          73 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSSVDSetDimensions_C",NULL));
     505          73 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSSVDGetDimensions_C",NULL));
     506          73 :   PetscFunctionReturn(PETSC_SUCCESS);
     507             : }
     508             : 
     509             : /*MC
     510             :    DSSVD - Dense Singular Value Decomposition.
     511             : 
     512             :    Level: beginner
     513             : 
     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.
     519             : 
     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().
     523             : 
     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.
     529             : 
     530             :    Used DS matrices:
     531             : +  DS_MAT_A - problem matrix
     532             : -  DS_MAT_T - upper bidiagonal matrix
     533             : 
     534             :    Implemented methods:
     535             : .  0 - Divide and Conquer (_bdsdc or _gesdd)
     536             : 
     537             : .seealso: DSCreate(), DSSetType(), DSType, DSSVDSetDimensions()
     538             : M*/
     539          73 : SLEPC_EXTERN PetscErrorCode DSCreate_SVD(DS ds)
     540             : {
     541          73 :   DS_SVD         *ctx;
     542             : 
     543          73 :   PetscFunctionBegin;
     544          73 :   PetscCall(PetscNew(&ctx));
     545          73 :   ds->data = (void*)ctx;
     546             : 
     547          73 :   ds->ops->allocate      = DSAllocate_SVD;
     548          73 :   ds->ops->view          = DSView_SVD;
     549          73 :   ds->ops->vectors       = DSVectors_SVD;
     550          73 :   ds->ops->solve[0]      = DSSolve_SVD_DC;
     551          73 :   ds->ops->sort          = DSSort_SVD;
     552             : #if !defined(PETSC_HAVE_MPIUNI)
     553          73 :   ds->ops->synchronize   = DSSynchronize_SVD;
     554             : #endif
     555          73 :   ds->ops->truncate      = DSTruncate_SVD;
     556          73 :   ds->ops->update        = DSUpdateExtraRow_SVD;
     557          73 :   ds->ops->destroy       = DSDestroy_SVD;
     558          73 :   ds->ops->matgetsize    = DSMatGetSize_SVD;
     559          73 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSSVDSetDimensions_C",DSSVDSetDimensions_SVD));
     560          73 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSSVDGetDimensions_C",DSSVDGetDimensions_SVD));
     561          73 :   PetscFunctionReturn(PETSC_SUCCESS);
     562             : }