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