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             :    This file implements a wrapper to the LAPACK SVD subroutines
      12             : */
      13             : 
      14             : #include <slepc/private/svdimpl.h>
      15             : #include <slepcblaslapack.h>
      16             : 
      17          25 : static PetscErrorCode SVDSetUp_LAPACK(SVD svd)
      18             : {
      19          25 :   PetscInt       M,N,P=0;
      20             : 
      21          25 :   PetscFunctionBegin;
      22          25 :   if (svd->nsv==0) svd->nsv = 1;
      23          25 :   PetscCall(MatGetSize(svd->A,&M,&N));
      24          25 :   if (!svd->isgeneralized) svd->ncv = N;
      25             :   else {
      26           7 :     PetscCall(MatGetSize(svd->OPb,&P,NULL));
      27           7 :     svd->ncv = PetscMin(M,PetscMin(N,P));
      28             :   }
      29          25 :   if (svd->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(svd,"Warning: parameter mpd ignored\n"));
      30          25 :   SVDCheckIgnored(svd,SVD_FEATURE_STOPPING);
      31          25 :   if (svd->max_it==PETSC_DETERMINE) svd->max_it = 1;
      32          25 :   svd->leftbasis = PETSC_TRUE;
      33          25 :   PetscCall(SVDAllocateSolution(svd,0));
      34          25 :   PetscCall(DSAllocate(svd->ds,PetscMax(N,PetscMax(M,P))));
      35          25 :   PetscFunctionReturn(PETSC_SUCCESS);
      36             : }
      37             : 
      38          16 : static PetscErrorCode SVDSolve_LAPACK(SVD svd)
      39             : {
      40          16 :   PetscInt          M,N,n,i,j,k,ld,lowu,lowv,highu,highv;
      41          16 :   Mat               A,Ar,mat;
      42          16 :   Vec               u,v;
      43          16 :   PetscScalar       *pU,*pV,*pu,*pv,*w;
      44             : 
      45          16 :   PetscFunctionBegin;
      46          16 :   PetscCall(DSGetLeadingDimension(svd->ds,&ld));
      47          16 :   PetscCall(MatCreateRedundantMatrix(svd->OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar));
      48          16 :   PetscCall(MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&mat));
      49          16 :   PetscCall(MatDestroy(&Ar));
      50          16 :   PetscCall(MatGetSize(mat,&M,&N));
      51          16 :   PetscCall(DSSetDimensions(svd->ds,M,0,0));
      52          16 :   PetscCall(DSSVDSetDimensions(svd->ds,N));
      53          16 :   PetscCall(DSGetMat(svd->ds,DS_MAT_A,&A));
      54          16 :   PetscCall(MatCopy(mat,A,SAME_NONZERO_PATTERN));
      55          16 :   PetscCall(DSRestoreMat(svd->ds,DS_MAT_A,&A));
      56          16 :   PetscCall(DSSetState(svd->ds,DS_STATE_RAW));
      57             : 
      58          16 :   n = PetscMin(M,N);
      59          16 :   PetscCall(PetscMalloc1(n,&w));
      60          16 :   PetscCall(DSSolve(svd->ds,w,NULL));
      61          16 :   PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
      62          16 :   PetscCall(DSSynchronize(svd->ds,w,NULL));
      63             : 
      64             :   /* copy singular vectors */
      65          16 :   PetscCall(DSGetArray(svd->ds,DS_MAT_U,&pU));
      66          16 :   PetscCall(DSGetArray(svd->ds,DS_MAT_V,&pV));
      67         454 :   for (i=0;i<n;i++) {
      68         438 :     if (svd->which == SVD_SMALLEST) k = n - i - 1;
      69             :     else k = i;
      70         438 :     svd->sigma[k] = PetscRealPart(w[i]);
      71         438 :     PetscCall(BVGetColumn(svd->U,k,&u));
      72         438 :     PetscCall(BVGetColumn(svd->V,k,&v));
      73         438 :     PetscCall(VecGetOwnershipRange(u,&lowu,&highu));
      74         438 :     PetscCall(VecGetOwnershipRange(v,&lowv,&highv));
      75         438 :     PetscCall(VecGetArray(u,&pu));
      76         438 :     PetscCall(VecGetArray(v,&pv));
      77         438 :     if (M>=N) {
      78       12492 :       for (j=lowu;j<highu;j++) pu[j-lowu] = pU[i*ld+j];
      79       11970 :       for (j=lowv;j<highv;j++) pv[j-lowv] = pV[i*ld+j];
      80             :     } else {
      81        1380 :       for (j=lowu;j<highu;j++) pu[j-lowu] = pV[i*ld+j];
      82        1260 :       for (j=lowv;j<highv;j++) pv[j-lowv] = pU[i*ld+j];
      83             :     }
      84         438 :     PetscCall(VecRestoreArray(u,&pu));
      85         438 :     PetscCall(VecRestoreArray(v,&pv));
      86         438 :     PetscCall(BVRestoreColumn(svd->U,k,&u));
      87         438 :     PetscCall(BVRestoreColumn(svd->V,k,&v));
      88             :   }
      89          16 :   PetscCall(DSRestoreArray(svd->ds,DS_MAT_U,&pU));
      90          16 :   PetscCall(DSRestoreArray(svd->ds,DS_MAT_V,&pV));
      91             : 
      92          16 :   svd->nconv  = n;
      93          16 :   svd->its    = 1;
      94          16 :   svd->reason = SVD_CONVERGED_TOL;
      95             : 
      96          16 :   PetscCall(MatDestroy(&mat));
      97          16 :   PetscCall(PetscFree(w));
      98          16 :   PetscFunctionReturn(PETSC_SUCCESS);
      99             : }
     100             : 
     101           7 : static PetscErrorCode SVDSolve_LAPACK_GSVD(SVD svd)
     102             : {
     103           7 :   PetscInt          nsv,m,n,p,i,j,mlocal,plocal,ld,lowx,lowu,lowv,highx;
     104           7 :   Mat               Ar,A,Ads,Br,B,Bds;
     105           7 :   Vec               uv,x;
     106           7 :   PetscScalar       *U,*V,*X,*px,*puv,*w;
     107             : 
     108           7 :   PetscFunctionBegin;
     109           7 :   PetscCall(DSGetLeadingDimension(svd->ds,&ld));
     110           7 :   PetscCall(MatCreateRedundantMatrix(svd->OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar));
     111           7 :   PetscCall(MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&A));
     112           7 :   PetscCall(MatDestroy(&Ar));
     113           7 :   PetscCall(MatCreateRedundantMatrix(svd->OPb,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Br));
     114           7 :   PetscCall(MatConvert(Br,MATSEQDENSE,MAT_INITIAL_MATRIX,&B));
     115           7 :   PetscCall(MatDestroy(&Br));
     116           7 :   PetscCall(MatGetSize(A,&m,&n));
     117           7 :   PetscCall(MatGetLocalSize(svd->OP,&mlocal,NULL));
     118           7 :   PetscCall(MatGetLocalSize(svd->OPb,&plocal,NULL));
     119           7 :   PetscCall(MatGetSize(B,&p,NULL));
     120           7 :   PetscCall(DSSetDimensions(svd->ds,m,0,0));
     121           7 :   PetscCall(DSGSVDSetDimensions(svd->ds,n,p));
     122           7 :   PetscCall(DSGetMat(svd->ds,DS_MAT_A,&Ads));
     123           7 :   PetscCall(MatCopy(A,Ads,SAME_NONZERO_PATTERN));
     124           7 :   PetscCall(DSRestoreMat(svd->ds,DS_MAT_A,&Ads));
     125           7 :   PetscCall(DSGetMat(svd->ds,DS_MAT_B,&Bds));
     126           7 :   PetscCall(MatCopy(B,Bds,SAME_NONZERO_PATTERN));
     127           7 :   PetscCall(DSRestoreMat(svd->ds,DS_MAT_B,&Bds));
     128           7 :   PetscCall(DSSetState(svd->ds,DS_STATE_RAW));
     129             : 
     130           7 :   nsv = PetscMin(n,PetscMin(p,m));
     131           7 :   PetscCall(PetscMalloc1(nsv,&w));
     132           7 :   PetscCall(DSSolve(svd->ds,w,NULL));
     133           7 :   PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
     134           7 :   PetscCall(DSSynchronize(svd->ds,w,NULL));
     135           7 :   PetscCall(DSGetDimensions(svd->ds,NULL,NULL,NULL,&nsv));
     136             : 
     137             :   /* copy singular vectors */
     138           7 :   PetscCall(MatGetOwnershipRange(svd->OP,&lowu,NULL));
     139           7 :   PetscCall(MatGetOwnershipRange(svd->OPb,&lowv,NULL));
     140           7 :   PetscCall(DSGetArray(svd->ds,DS_MAT_X,&X));
     141           7 :   PetscCall(DSGetArray(svd->ds,DS_MAT_U,&U));
     142           7 :   PetscCall(DSGetArray(svd->ds,DS_MAT_V,&V));
     143          98 :   for (j=0;j<nsv;j++) {
     144          91 :     svd->sigma[j] = PetscRealPart(w[j]);
     145          91 :     PetscCall(BVGetColumn(svd->V,j,&x));
     146          91 :     PetscCall(VecGetOwnershipRange(x,&lowx,&highx));
     147          91 :     PetscCall(VecGetArrayWrite(x,&px));
     148        1776 :     for (i=lowx;i<highx;i++) px[i-lowx] = X[i+j*ld];
     149          91 :     PetscCall(VecRestoreArrayWrite(x,&px));
     150          91 :     PetscCall(BVRestoreColumn(svd->V,j,&x));
     151          91 :     PetscCall(BVGetColumn(svd->U,j,&uv));
     152          91 :     PetscCall(VecGetArrayWrite(uv,&puv));
     153        1701 :     for (i=0;i<mlocal;i++) puv[i] = U[i+lowu+j*ld];
     154        1877 :     for (i=0;i<plocal;i++) puv[i+mlocal] = V[i+lowv+j*ld];
     155          91 :     PetscCall(VecRestoreArrayWrite(uv,&puv));
     156          91 :     PetscCall(BVRestoreColumn(svd->U,j,&uv));
     157             :   }
     158           7 :   PetscCall(DSRestoreArray(svd->ds,DS_MAT_X,&X));
     159           7 :   PetscCall(DSRestoreArray(svd->ds,DS_MAT_U,&U));
     160           7 :   PetscCall(DSRestoreArray(svd->ds,DS_MAT_V,&V));
     161             : 
     162           7 :   svd->nconv  = nsv;
     163           7 :   svd->its    = 1;
     164           7 :   svd->reason = SVD_CONVERGED_TOL;
     165             : 
     166           7 :   PetscCall(MatDestroy(&A));
     167           7 :   PetscCall(MatDestroy(&B));
     168           7 :   PetscCall(PetscFree(w));
     169           7 :   PetscFunctionReturn(PETSC_SUCCESS);
     170             : }
     171             : 
     172           2 : static PetscErrorCode SVDSolve_LAPACK_HSVD(SVD svd)
     173             : {
     174           2 :   PetscInt          M,N,n,i,j,k,ld,lowu,lowv,highu,highv;
     175           2 :   Mat               A,Ar,mat,D;
     176           2 :   Vec               u,v,vomega;
     177           2 :   PetscScalar       *pU,*pV,*pu,*pv,*w;
     178           2 :   PetscReal         *pD;
     179             : 
     180           2 :   PetscFunctionBegin;
     181           2 :   PetscCall(DSGetLeadingDimension(svd->ds,&ld));
     182           2 :   PetscCall(MatCreateRedundantMatrix(svd->OP,0,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&Ar));
     183           2 :   PetscCall(MatConvert(Ar,MATSEQDENSE,MAT_INITIAL_MATRIX,&mat));
     184           2 :   PetscCall(MatDestroy(&Ar));
     185           2 :   PetscCall(MatGetSize(mat,&M,&N));
     186           2 :   PetscCall(DSSetDimensions(svd->ds,M,0,0));
     187           2 :   PetscCall(DSHSVDSetDimensions(svd->ds,N));
     188           2 :   PetscCall(DSGetMat(svd->ds,DS_MAT_A,&A));
     189           2 :   PetscCall(MatCopy(mat,A,SAME_NONZERO_PATTERN));
     190           2 :   PetscCall(DSRestoreMat(svd->ds,DS_MAT_A,&A));
     191           2 :   PetscCall(DSGetMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
     192           2 :   PetscCall(VecCopy(svd->omega,vomega));
     193           2 :   PetscCall(DSRestoreMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
     194           2 :   PetscCall(DSSetState(svd->ds,DS_STATE_RAW));
     195             : 
     196           2 :   n = PetscMin(M,N);
     197           2 :   PetscCall(PetscMalloc1(n,&w));
     198           2 :   PetscCall(DSSolve(svd->ds,w,NULL));
     199           2 :   PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
     200           2 :   PetscCall(DSSynchronize(svd->ds,w,NULL));
     201             : 
     202             :   /* copy singular vectors */
     203           2 :   PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&pD));
     204           2 :   PetscCall(DSGetArray(svd->ds,DS_MAT_U,&pU));
     205           2 :   PetscCall(DSGetArray(svd->ds,DS_MAT_V,&pV));
     206         522 :   for (i=0;i<n;i++) {
     207         520 :     if (svd->which == SVD_SMALLEST) k = n - i - 1;
     208             :     else k = i;
     209         520 :     svd->sigma[k] = PetscRealPart(w[i]);
     210         520 :     svd->sign[k]  = pD[i];
     211         520 :     PetscCall(BVGetColumn(svd->U,k,&u));
     212         520 :     PetscCall(BVGetColumn(svd->V,k,&v));
     213         520 :     PetscCall(VecGetOwnershipRange(u,&lowu,&highu));
     214         520 :     PetscCall(VecGetOwnershipRange(v,&lowv,&highv));
     215         520 :     PetscCall(VecGetArray(u,&pu));
     216         520 :     PetscCall(VecGetArray(v,&pv));
     217         520 :     if (M>=N) {
     218      371080 :       for (j=lowu;j<highu;j++) pu[j-lowu] = pU[i*ld+j];
     219      142920 :       for (j=lowv;j<highv;j++) pv[j-lowv] = pV[i*ld+j];
     220             :     } else {
     221           0 :       for (j=lowu;j<highu;j++) pu[j-lowu] = pV[i*ld+j];
     222           0 :       for (j=lowv;j<highv;j++) pv[j-lowv] = pU[i*ld+j];
     223             :     }
     224         520 :     PetscCall(VecRestoreArray(u,&pu));
     225         520 :     PetscCall(VecRestoreArray(v,&pv));
     226         520 :     PetscCall(BVRestoreColumn(svd->U,k,&u));
     227         520 :     PetscCall(BVRestoreColumn(svd->V,k,&v));
     228             :   }
     229           2 :   PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&pD));
     230           2 :   PetscCall(DSRestoreArray(svd->ds,DS_MAT_U,&pU));
     231           2 :   PetscCall(DSRestoreArray(svd->ds,DS_MAT_V,&pV));
     232             : 
     233           2 :   svd->nconv  = n;
     234           2 :   svd->its    = 1;
     235           2 :   svd->reason = SVD_CONVERGED_TOL;
     236             : 
     237           2 :   PetscCall(MatDestroy(&mat));
     238           2 :   PetscCall(PetscFree(w));
     239           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     240             : }
     241             : 
     242          43 : static PetscErrorCode SVDSetDSType_LAPACK(SVD svd)
     243             : {
     244          43 :   PetscFunctionBegin;
     245          69 :   PetscCall(DSSetType(svd->ds,svd->OPb?DSGSVD:svd->omega?DSHSVD:DSSVD));
     246          43 :   PetscFunctionReturn(PETSC_SUCCESS);
     247             : }
     248             : 
     249          19 : SLEPC_EXTERN PetscErrorCode SVDCreate_LAPACK(SVD svd)
     250             : {
     251          19 :   PetscFunctionBegin;
     252          19 :   svd->ops->setup     = SVDSetUp_LAPACK;
     253          19 :   svd->ops->solve     = SVDSolve_LAPACK;
     254          19 :   svd->ops->solveg    = SVDSolve_LAPACK_GSVD;
     255          19 :   svd->ops->solveh    = SVDSolve_LAPACK_HSVD;
     256          19 :   svd->ops->setdstype = SVDSetDSType_LAPACK;
     257          19 :   PetscFunctionReturn(PETSC_SUCCESS);
     258             : }