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             :    Common subroutines for all Krylov-type PEP solvers
      12             : */
      13             : 
      14             : #include <slepc/private/pepimpl.h>    /*I "slepcpep.h" I*/
      15             : #include <slepcblaslapack.h>
      16             : #include "pepkrylov.h"
      17             : 
      18         103 : PetscErrorCode PEPExtractVectors_TOAR(PEP pep)
      19             : {
      20         103 :   PetscInt          i,j,nq,deg=pep->nmat-1,lds,idxcpy=0,ldds,k,ld;
      21         103 :   PetscScalar       *X,*er,*ei,*SS,*vals,*ivals,sone=1.0,szero=0.0,*yi,*yr,*tr,*ti,alpha,*pS0;
      22         103 :   const PetscScalar *S;
      23         103 :   PetscBLASInt      k_,nq_,lds_,one=1,ldds_,cols,info,zero=0;
      24         103 :   PetscBool         flg;
      25         103 :   PetscReal         norm,max,t,factor=1.0,done=1.0;
      26         103 :   Vec               xr,xi,w[4];
      27         103 :   PEP_TOAR          *ctx = (PEP_TOAR*)pep->data;
      28         103 :   Mat               S0,MS;
      29             : 
      30         103 :   PetscFunctionBegin;
      31         103 :   PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
      32         103 :   PetscCall(MatDenseGetArrayRead(MS,&S));
      33         103 :   PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
      34         103 :   PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
      35         103 :   k = pep->nconv;
      36         103 :   if (k==0) PetscFunctionReturn(PETSC_SUCCESS);
      37         103 :   lds = deg*ld;
      38         103 :   PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
      39         103 :   PetscCall(PetscCalloc5(k,&er,k,&ei,nq*k,&SS,pep->nmat,&vals,pep->nmat,&ivals));
      40         103 :   PetscCall(STGetTransform(pep->st,&flg));
      41         103 :   if (flg) factor = pep->sfactor;
      42         700 :   for (i=0;i<k;i++) {
      43         597 :     er[i] = factor*pep->eigr[i];
      44         597 :     ei[i] = factor*pep->eigi[i];
      45             :   }
      46         103 :   PetscCall(STBackTransform(pep->st,k,er,ei));
      47             : 
      48         103 :   PetscCall(DSVectors(pep->ds,DS_MAT_X,NULL,NULL));
      49         103 :   PetscCall(DSGetArray(pep->ds,DS_MAT_X,&X));
      50             : 
      51         103 :   PetscCall(PetscBLASIntCast(k,&k_));
      52         103 :   PetscCall(PetscBLASIntCast(nq,&nq_));
      53         103 :   PetscCall(PetscBLASIntCast(lds,&lds_));
      54         103 :   PetscCall(PetscBLASIntCast(ldds,&ldds_));
      55             : 
      56         103 :   if (pep->extract==PEP_EXTRACT_NONE || pep->refine==PEP_REFINE_MULTIPLE) {
      57          36 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&nq_,&k_,&k_,&sone,S,&lds_,X,&ldds_,&szero,SS,&nq_));
      58             :   } else {
      59          67 :     switch (pep->extract) {
      60             :     case PEP_EXTRACT_NONE:
      61             :       break;
      62             :     case PEP_EXTRACT_NORM:
      63         341 :       for (i=0;i<k;i++) {
      64         279 :         PetscCall(PEPEvaluateBasis(pep,er[i],ei[i],vals,ivals));
      65             :         max = 1.0;
      66         577 :         for (j=1;j<deg;j++) {
      67         298 :           norm = SlepcAbsEigenvalue(vals[j],ivals[j]);
      68         298 :           if (max < norm) { max = norm; idxcpy = j; }
      69             :         }
      70         279 :         PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
      71             : #if !defined(PETSC_USE_COMPLEX)
      72         279 :         if (PetscRealPart(ei[i])!=0.0) {
      73          44 :           i++;
      74         279 :           PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
      75             :         }
      76             : #endif
      77             :       }
      78             :       break;
      79           3 :     case PEP_EXTRACT_RESIDUAL:
      80           3 :       PetscCall(VecDuplicate(pep->work[0],&xr));
      81           3 :       PetscCall(VecDuplicate(pep->work[0],&w[0]));
      82           3 :       PetscCall(VecDuplicate(pep->work[0],&w[1]));
      83             : #if !defined(PETSC_USE_COMPLEX)
      84           3 :       PetscCall(VecDuplicate(pep->work[0],&w[2]));
      85           3 :       PetscCall(VecDuplicate(pep->work[0],&w[3]));
      86           3 :       PetscCall(VecDuplicate(pep->work[0],&xi));
      87             : #else
      88             :       xi = NULL;
      89             : #endif
      90          23 :       for (i=0;i<k;i++) {
      91             :         max = 0.0;
      92          63 :         for (j=0;j<deg;j++) {
      93          43 :           PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+j*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
      94          43 :           PetscCall(BVMultVec(pep->V,1.0,0.0,xr,SS+i*nq));
      95             : #if !defined(PETSC_USE_COMPLEX)
      96          43 :           PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+j*ld,&lds_,X+(i+1)*ldds,&one,&szero,SS+i*nq,&one));
      97          43 :           PetscCall(BVMultVec(pep->V,1.0,0.0,xi,SS+i*nq));
      98             : #endif
      99          43 :           PetscCall(PEPComputeResidualNorm_Private(pep,er[i],ei[i],xr,xi,w,&norm));
     100          43 :           if (norm>max) { max = norm; idxcpy=j; }
     101             :         }
     102          20 :         PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
     103             : #if !defined(PETSC_USE_COMPLEX)
     104          20 :         if (PetscRealPart(ei[i])!=0.0) {
     105           3 :           i++;
     106          20 :           PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
     107             :         }
     108             : #endif
     109             :       }
     110           3 :       PetscCall(VecDestroy(&xr));
     111           3 :       PetscCall(VecDestroy(&w[0]));
     112           3 :       PetscCall(VecDestroy(&w[1]));
     113             : #if !defined(PETSC_USE_COMPLEX)
     114           3 :       PetscCall(VecDestroy(&w[2]));
     115           3 :       PetscCall(VecDestroy(&w[3]));
     116           3 :       PetscCall(VecDestroy(&xi));
     117             : #endif
     118             :       break;
     119           2 :     case PEP_EXTRACT_STRUCTURED:
     120           2 :       PetscCall(PetscMalloc2(k,&tr,k,&ti));
     121          20 :       for (i=0;i<k;i++) {
     122          18 :         t = 0.0;
     123          18 :         PetscCall(PEPEvaluateBasis(pep,er[i],ei[i],vals,ivals));
     124          18 :         yr = X+i*ldds; yi = NULL;
     125             : #if !defined(PETSC_USE_COMPLEX)
     126          18 :         if (ei[i]!=0.0) { yr = tr; yi = ti; }
     127             : #endif
     128          54 :         for (j=0;j<deg;j++) {
     129          36 :           alpha = PetscConj(vals[j]);
     130             : #if !defined(PETSC_USE_COMPLEX)
     131          36 :           if (ei[i]!=0.0) {
     132           0 :             PetscCall(PetscArrayzero(tr,k));
     133           0 :             PetscCallBLAS("BLASaxpy",BLASaxpy_(&k_,&vals[j],X+i*ldds,&one,tr,&one));
     134           0 :             PetscCallBLAS("BLASaxpy",BLASaxpy_(&k_,&ivals[j],X+(i+1)*ldds,&one,tr,&one));
     135           0 :             PetscCall(PetscArrayzero(ti,k));
     136           0 :             PetscCallBLAS("BLASaxpy",BLASaxpy_(&k_,&vals[j],X+(i+1)*ldds,&one,ti,&one));
     137           0 :             alpha = -ivals[j];
     138           0 :             PetscCallBLAS("BLASaxpy",BLASaxpy_(&k_,&alpha,X+i*ldds,&one,ti,&one));
     139           0 :             alpha = 1.0;
     140             :           }
     141             : #endif
     142          36 :           PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&alpha,S+j*ld,&lds_,yr,&one,&sone,SS+i*nq,&one));
     143          36 :           t += SlepcAbsEigenvalue(vals[j],ivals[j])*SlepcAbsEigenvalue(vals[j],ivals[j]);
     144          36 :           if (yi) {
     145          36 :             PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&alpha,S+j*ld,&lds_,yi,&one,&sone,SS+(i+1)*nq,&one));
     146             :           }
     147             :         }
     148          18 :         cols = yi? 2: 1;
     149          18 :         PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&t,&done,&nq_,&cols,SS+i*nq,&nq_,&info));
     150          18 :         SlepcCheckLapackInfo("lascl",info);
     151          18 :         if (yi) i++;
     152             :       }
     153           2 :       PetscCall(PetscFree2(tr,ti));
     154             :       break;
     155             :     }
     156         103 :   }
     157             : 
     158             :   /* update vectors V = V*S */
     159         103 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nq,k,NULL,&S0));
     160         103 :   PetscCall(MatDenseGetArrayWrite(S0,&pS0));
     161         700 :   for (i=0;i<k;i++) PetscCall(PetscArraycpy(pS0+i*nq,SS+i*nq,nq));
     162         103 :   PetscCall(MatDenseRestoreArrayWrite(S0,&pS0));
     163         103 :   PetscCall(BVMultInPlace(pep->V,S0,0,k));
     164         103 :   PetscCall(MatDestroy(&S0));
     165         103 :   PetscCall(PetscFree5(er,ei,SS,vals,ivals));
     166         103 :   PetscCall(MatDenseRestoreArrayRead(MS,&S));
     167         103 :   PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
     168         103 :   PetscFunctionReturn(PETSC_SUCCESS);
     169             : }
     170             : 
     171             : /*
     172             :    PEPKrylovConvergence - This is the analogue to EPSKrylovConvergence, but
     173             :    for polynomial Krylov methods.
     174             : 
     175             :    Differences:
     176             :    - Always non-symmetric
     177             :    - Does not check for STSHIFT
     178             :    - No correction factor
     179             :    - No support for true residual
     180             : */
     181        1208 : PetscErrorCode PEPKrylovConvergence(PEP pep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscInt *kout)
     182             : {
     183        1208 :   PetscInt       k,newk,marker,inside;
     184        1208 :   PetscScalar    re,im;
     185        1208 :   PetscReal      resnorm;
     186        1208 :   PetscBool      istrivial;
     187             : 
     188        1208 :   PetscFunctionBegin;
     189        1208 :   PetscCall(RGIsTrivial(pep->rg,&istrivial));
     190        1208 :   marker = -1;
     191        1208 :   if (pep->trackall) getall = PETSC_TRUE;
     192        2144 :   for (k=kini;k<kini+nits;k++) {
     193             :     /* eigenvalue */
     194        2142 :     re = pep->eigr[k];
     195        2142 :     im = pep->eigi[k];
     196        2142 :     if (PetscUnlikely(!istrivial)) {
     197         395 :       PetscCall(STBackTransform(pep->st,1,&re,&im));
     198         395 :       PetscCall(RGCheckInside(pep->rg,1,&re,&im,&inside));
     199         395 :       if (marker==-1 && inside<0) marker = k;
     200         395 :       re = pep->eigr[k];
     201         395 :       im = pep->eigi[k];
     202             :     }
     203        2142 :     newk = k;
     204        2142 :     PetscCall(DSVectors(pep->ds,DS_MAT_X,&newk,&resnorm));
     205        2142 :     resnorm *= beta;
     206             :     /* error estimate */
     207        2142 :     PetscCall((*pep->converged)(pep,re,im,resnorm,&pep->errest[k],pep->convergedctx));
     208        2142 :     if (marker==-1 && pep->errest[k] >= pep->tol) marker = k;
     209        2142 :     if (PetscUnlikely(newk==k+1)) {
     210         281 :       pep->errest[k+1] = pep->errest[k];
     211         281 :       k++;
     212             :     }
     213        2142 :     if (marker!=-1 && !getall) break;
     214             :   }
     215        1208 :   if (marker!=-1) k = marker;
     216        1208 :   *kout = k;
     217        1208 :   PetscFunctionReturn(PETSC_SUCCESS);
     218             : }