Subversion Repositories slepc-dev

Rev

Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
1376 slepc 1 
/*
6 dsic.upv.es!jroman 2 
       This file implements a wrapper to the LAPACK eigenvalue subroutines.
679 dsic.upv.es!antodo 3 
       Generalized problems are transformed to standard ones only if necessary.
1376 slepc 4 
 
5 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1672 slepc 6 
   SLEPc - Scalable Library for Eigenvalue Problem Computations
2116 eromero 7 
   Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
1376 slepc 8 
 
1672 slepc 9 
   This file is part of SLEPc.
10 
 
11 
   SLEPc is free software: you can redistribute it and/or modify it under  the
12 
   terms of version 3 of the GNU Lesser General Public License as published by
13 
   the Free Software Foundation.
14 
 
15 
   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY
16 
   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS
17 
   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for
18 
   more details.
19 
 
20 
   You  should have received a copy of the GNU Lesser General  Public  License
21 
   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
1376 slepc 22 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
6 dsic.upv.es!jroman 23 
*/
1376 slepc 24 
 
2283 jroman 25 
#include <private/epsimpl.h>     /*I "slepceps.h" I*/
26 
#include <slepcblaslapack.h>
6 dsic.upv.es!jroman 27 
 
679 dsic.upv.es!antodo 28 
typedef struct {
29 
  Mat OP,A,B;
30 
} EPS_LAPACK;
31 
 
6 dsic.upv.es!jroman 32 
#undef __FUNCT__  
33 
#define __FUNCT__ "EPSSetUp_LAPACK"
476 dsic.upv.es!antodo 34 
PetscErrorCode EPSSetUp_LAPACK(EPS eps)
6 dsic.upv.es!jroman 35 
{
1096 slepc 36 
  PetscErrorCode ierr,ierra,ierrb;
476 dsic.upv.es!antodo 37 
  EPS_LAPACK     *la = (EPS_LAPACK *)eps->data;
2216 jroman 38 
  PetscBool      flg;
679 dsic.upv.es!antodo 39 
  Mat            A,B;
40 
  PetscScalar    shift;
1245 slepc 41 
  KSP            ksp;
42 
  PC             pc;
679 dsic.upv.es!antodo 43 
 
6 dsic.upv.es!jroman 44 
  PetscFunctionBegin;
1928 jroman 45 
  eps->ncv = eps->n;
1579 slepc 46 
  if (eps->mpd) PetscInfo(eps,"Warning: parameter mpd ignored\n");
260 dsic.upv.es!antodo 47 
 
1942 jroman 48 
  if (!eps->which) eps->which = EPS_LARGEST_MAGNITUDE;
1940 jroman 49 
  if (eps->balance!=EPS_BALANCE_NONE)
2214 jroman 50 
    SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Balancing not supported in Lapack solvers");
1819 jroman 51 
 
2306 jroman 52 
  ierr = MatDestroy(&la->OP);CHKERRQ(ierr);
53 
  ierr = MatDestroy(&la->A);CHKERRQ(ierr);
54 
  ierr = MatDestroy(&la->B);CHKERRQ(ierr);
679 dsic.upv.es!antodo 55 
 
56 
  ierr = PetscTypeCompare((PetscObject)eps->OP,STSHIFT,&flg);CHKERRQ(ierr);
1079 slepc 57 
  ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
58 
 
1096 slepc 59 
  if (flg) {
679 dsic.upv.es!antodo 60 
    la->OP = PETSC_NULL;
1245 slepc 61 
    PetscPushErrorHandler(PetscIgnoreErrorHandler,PETSC_NULL);
1096 slepc 62 
    ierra = SlepcMatConvertSeqDense(A,&la->A);CHKERRQ(ierr);
679 dsic.upv.es!antodo 63 
    if (eps->isgeneralized) {
1096 slepc 64 
      ierrb = SlepcMatConvertSeqDense(B,&la->B);CHKERRQ(ierr);
65 
    } else {
66 
      ierrb = 0;
67 
      la->B = PETSC_NULL;
679 dsic.upv.es!antodo 68 
    }
1096 slepc 69 
    PetscPopErrorHandler();
70 
    if (ierra == 0 && ierrb == 0) {
71 
      ierr = STGetShift(eps->OP,&shift);CHKERRQ(ierr);
72 
      if (shift != 0.0) {
2316 jroman 73 
        ierr = MatShift(la->A,shift);CHKERRQ(ierr);
1096 slepc 74 
      }
1245 slepc 75 
      /* use dummy pc and ksp to avoid problems when B is not positive definite */
76 
      ierr = STGetKSP(eps->OP,&ksp);CHKERRQ(ierr);
77 
      ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
78 
      ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
79 
      ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
1596 slepc 80 
      ierr = EPSAllocateSolution(eps);CHKERRQ(ierr);
1096 slepc 81 
      PetscFunctionReturn(0);
679 dsic.upv.es!antodo 82 
    }
6 dsic.upv.es!jroman 83 
  }
1096 slepc 84 
  PetscInfo(eps,"Using slow explicit operator\n");
85 
  la->A = PETSC_NULL;
86 
  la->B = PETSC_NULL;
87 
  ierr = STComputeExplicitOperator(eps->OP,&la->OP);CHKERRQ(ierr);
88 
  ierr = PetscTypeCompare((PetscObject)la->OP,MATSEQDENSE,&flg);CHKERRQ(ierr);
89 
  if (!flg) {
90 
    ierr = SlepcMatConvertSeqDense(la->OP,&la->OP);CHKERRQ(ierr);
91 
  }
1560 slepc 92 
  if (eps->extraction) {
93 
     ierr = PetscInfo(eps,"Warning: extraction type ignored\n");CHKERRQ(ierr);
1426 slepc 94 
  }
1596 slepc 95 
  ierr = EPSAllocateSolution(eps);CHKERRQ(ierr);
6 dsic.upv.es!jroman 96 
  PetscFunctionReturn(0);
97 
}
98 
 
99 
#undef __FUNCT__  
100 
#define __FUNCT__ "EPSSolve_LAPACK"
476 dsic.upv.es!antodo 101 
PetscErrorCode EPSSolve_LAPACK(EPS eps)
6 dsic.upv.es!jroman 102 
{
476 dsic.upv.es!antodo 103 
  PetscErrorCode ierr;
1928 jroman 104 
  PetscInt       n=eps->n,i,low,high;
982 slepc 105 
  PetscMPIInt    size;
780 dsic.upv.es!jroman 106 
  PetscScalar    *array,*arrayb,*pV,*pW;
679 dsic.upv.es!antodo 107 
  PetscReal      *w;
476 dsic.upv.es!antodo 108 
  EPS_LAPACK     *la = (EPS_LAPACK *)eps->data;
1422 slepc 109 
  MPI_Comm       comm = ((PetscObject)eps)->comm;
6 dsic.upv.es!jroman 110 
 
111 
  PetscFunctionBegin;
112 
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
393 dsic.upv.es!antodo 113 
  if (size == 1) {
114 
    ierr = VecGetArray(eps->V[0],&pV);CHKERRQ(ierr);
115 
  } else {
679 dsic.upv.es!antodo 116 
    ierr = PetscMalloc(sizeof(PetscScalar)*n*n,&pV);CHKERRQ(ierr);
393 dsic.upv.es!antodo 117 
  }
1947 jroman 118 
  if (eps->leftvecs) {
780 dsic.upv.es!jroman 119 
    if (size == 1) {
120 
      ierr = VecGetArray(eps->W[0],&pW);CHKERRQ(ierr);
121 
    } else {
122 
      ierr = PetscMalloc(sizeof(PetscScalar)*n*n,&pW);CHKERRQ(ierr);
123 
    }
124 
  } else pW = PETSC_NULL;
73 dsic.upv.es!antodo 125 
 
679 dsic.upv.es!antodo 126 
  if (la->OP) {
127 
    ierr = MatGetArray(la->OP,&array);CHKERRQ(ierr);
780 dsic.upv.es!jroman 128 
    ierr = EPSDenseNHEP(n,array,eps->eigr,eps->eigi,pV,pW);CHKERRQ(ierr);  
679 dsic.upv.es!antodo 129 
    ierr = MatRestoreArray(la->OP,&array);CHKERRQ(ierr);
130 
  } else if (eps->ishermitian) {
131 
#if defined(PETSC_USE_COMPLEX)
132 
    ierr = PetscMalloc(n*sizeof(PetscReal),&w);CHKERRQ(ierr);
133 
#else
134 
    w = eps->eigr;
135 
#endif
136 
    ierr = MatGetArray(la->A,&array);CHKERRQ(ierr);
137 
    if (!eps->isgeneralized) {
1175 slepc 138 
      ierr = EPSDenseHEP(n,array,n,w,pV);CHKERRQ(ierr);
679 dsic.upv.es!antodo 139 
    } else {
140 
      ierr = MatGetArray(la->B,&arrayb);CHKERRQ(ierr);
141 
      ierr = EPSDenseGHEP(n,array,arrayb,w,pV);CHKERRQ(ierr);  
142 
      ierr = MatRestoreArray(la->B,&arrayb);CHKERRQ(ierr);
143 
    }
144 
    ierr = MatRestoreArray(la->A,&array);CHKERRQ(ierr);
145 
#if defined(PETSC_USE_COMPLEX)
146 
    for (i=0;i<n;i++) {
147 
      eps->eigr[i] = w[i];
148 
    }
149 
    ierr = PetscFree(w);CHKERRQ(ierr);
150 
#endif    
151 
  } else {
152 
    ierr = MatGetArray(la->A,&array);CHKERRQ(ierr);
153 
    if (!eps->isgeneralized) {
780 dsic.upv.es!jroman 154 
      ierr = EPSDenseNHEP(n,array,eps->eigr,eps->eigi,pV,pW);CHKERRQ(ierr);
679 dsic.upv.es!antodo 155 
    } else {
156 
      ierr = MatGetArray(la->B,&arrayb);CHKERRQ(ierr);
780 dsic.upv.es!jroman 157 
      ierr = EPSDenseGNHEP(n,array,arrayb,eps->eigr,eps->eigi,pV,pW);CHKERRQ(ierr);  
679 dsic.upv.es!antodo 158 
      ierr = MatRestoreArray(la->B,&arrayb);CHKERRQ(ierr);
159 
    }
160 
    ierr = MatRestoreArray(la->A,&array);CHKERRQ(ierr);
161 
  }
6 dsic.upv.es!jroman 162 
 
393 dsic.upv.es!antodo 163 
  if (size == 1) {
164 
    ierr = VecRestoreArray(eps->V[0],&pV);CHKERRQ(ierr);
165 
  } else {
166 
    for (i=0; i<eps->ncv; i++) {
167 
      ierr = VecGetOwnershipRange(eps->V[i], &low, &high);CHKERRQ(ierr);
168 
      ierr = VecGetArray(eps->V[i], &array);CHKERRQ(ierr);
169 
      ierr = PetscMemcpy(array, pV+i*n+low, (high-low)*sizeof(PetscScalar));
170 
      ierr = VecRestoreArray(eps->V[i], &array);CHKERRQ(ierr);
73 dsic.upv.es!antodo 171 
    }
393 dsic.upv.es!antodo 172 
    ierr = PetscFree(pV);CHKERRQ(ierr);
6 dsic.upv.es!jroman 173 
  }
780 dsic.upv.es!jroman 174 
  if (pW) {
175 
    if (size == 1) {
176 
      ierr = VecRestoreArray(eps->W[0],&pW);CHKERRQ(ierr);
177 
    } else {
178 
      for (i=0; i<eps->ncv; i++) {
179 
        ierr = VecGetOwnershipRange(eps->W[i], &low, &high);CHKERRQ(ierr);
180 
        ierr = VecGetArray(eps->W[i], &array);CHKERRQ(ierr);
181 
        ierr = PetscMemcpy(array, pW+i*n+low, (high-low)*sizeof(PetscScalar));
182 
        ierr = VecRestoreArray(eps->W[i], &array);CHKERRQ(ierr);
183 
      }
184 
      ierr = PetscFree(pW);CHKERRQ(ierr);
185 
    }
186 
  }
6 dsic.upv.es!jroman 187 
  eps->nconv = eps->ncv;
146 dsic.upv.es!antodo 188 
  eps->its   = 1;  
189 
  eps->reason = EPS_CONVERGED_TOL;
6 dsic.upv.es!jroman 190 
  PetscFunctionReturn(0);
191 
}
192 
 
193 
#undef __FUNCT__  
194 
#define __FUNCT__ "EPSDestroy_LAPACK"
476 dsic.upv.es!antodo 195 
PetscErrorCode EPSDestroy_LAPACK(EPS eps)
6 dsic.upv.es!jroman 196 
{
476 dsic.upv.es!antodo 197 
  PetscErrorCode ierr;
198 
  EPS_LAPACK     *la = (EPS_LAPACK *)eps->data;
6 dsic.upv.es!jroman 199 
 
200 
  PetscFunctionBegin;
2213 jroman 201 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2306 jroman 202 
  ierr = MatDestroy(&la->OP);CHKERRQ(ierr);
203 
  ierr = MatDestroy(&la->A);CHKERRQ(ierr);
204 
  ierr = MatDestroy(&la->B);CHKERRQ(ierr);
1040 slepc 205 
  ierr = PetscFree(eps->data);CHKERRQ(ierr);
1596 slepc 206 
  ierr = EPSFreeSolution(eps);CHKERRQ(ierr);
6 dsic.upv.es!jroman 207 
  PetscFunctionReturn(0);
208 
}
209 
 
210 
EXTERN_C_BEGIN
211 
#undef __FUNCT__  
212 
#define __FUNCT__ "EPSCreate_LAPACK"
476 dsic.upv.es!antodo 213 
PetscErrorCode EPSCreate_LAPACK(EPS eps)
6 dsic.upv.es!jroman 214 
{
476 dsic.upv.es!antodo 215 
  PetscErrorCode ierr;
6 dsic.upv.es!jroman 216 
 
217 
  PetscFunctionBegin;
2329 jroman 218 
  ierr = PetscNewLog(eps,EPS_LAPACK,&eps->data);CHKERRQ(ierr);
6 dsic.upv.es!jroman 219 
  eps->ops->solve                = EPSSolve_LAPACK;
503 dsic.upv.es!antodo 220 
  eps->ops->setup                = EPSSetUp_LAPACK;
6 dsic.upv.es!jroman 221 
  eps->ops->destroy              = EPSDestroy_LAPACK;
177 dsic.upv.es!antodo 222 
  eps->ops->backtransform        = EPSBackTransform_Default;
503 dsic.upv.es!antodo 223 
  eps->ops->computevectors       = EPSComputeVectors_Default;
6 dsic.upv.es!jroman 224 
  PetscFunctionReturn(0);
225 
}
226 
EXTERN_C_END