Subversion Repositories slepc-dev

Rev

Go to most recent revision | 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
7 
   Copyright (c) 2002-2009, 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 
 
1521 slepc 25 
#include "private/epsimpl.h"
6 dsic.upv.es!jroman 26 
#include "slepcblaslapack.h"
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;
1096 slepc 38 
  PetscTruth     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 
 
1819 jroman 48 
  if (eps->balance!=EPSBALANCE_NONE)
49 
    SETERRQ(PETSC_ERR_SUP,"Balancing not supported in Lapack solvers");
50 
 
679 dsic.upv.es!antodo 51 
  if (la->OP) { ierr = MatDestroy(la->OP);CHKERRQ(ierr); }
52 
  if (la->A) { ierr = MatDestroy(la->A);CHKERRQ(ierr); }
53 
  if (la->B) { ierr = MatDestroy(la->B);CHKERRQ(ierr); }
54 
 
55 
  ierr = PetscTypeCompare((PetscObject)eps->OP,STSHIFT,&flg);CHKERRQ(ierr);
1079 slepc 56 
  ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
57 
 
1096 slepc 58 
  if (flg) {
679 dsic.upv.es!antodo 59 
    la->OP = PETSC_NULL;
1245 slepc 60 
    PetscPushErrorHandler(PetscIgnoreErrorHandler,PETSC_NULL);
1096 slepc 61 
    ierra = SlepcMatConvertSeqDense(A,&la->A);CHKERRQ(ierr);
679 dsic.upv.es!antodo 62 
    if (eps->isgeneralized) {
1096 slepc 63 
      ierrb = SlepcMatConvertSeqDense(B,&la->B);CHKERRQ(ierr);
64 
    } else {
65 
      ierrb = 0;
66 
      la->B = PETSC_NULL;
679 dsic.upv.es!antodo 67 
    }
1096 slepc 68 
    PetscPopErrorHandler();
69 
    if (ierra == 0 && ierrb == 0) {
70 
      ierr = STGetShift(eps->OP,&shift);CHKERRQ(ierr);
71 
      if (shift != 0.0) {
72 
        ierr = MatShift(la->A,shift);CHKERRQ(ierr);
73 
      }
1245 slepc 74 
      /* use dummy pc and ksp to avoid problems when B is not positive definite */
75 
      ierr = STGetKSP(eps->OP,&ksp);CHKERRQ(ierr);
76 
      ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
77 
      ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
78 
      ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
1596 slepc 79 
      ierr = EPSAllocateSolution(eps);CHKERRQ(ierr);
1096 slepc 80 
      PetscFunctionReturn(0);
679 dsic.upv.es!antodo 81 
    }
6 dsic.upv.es!jroman 82 
  }
1096 slepc 83 
  PetscInfo(eps,"Using slow explicit operator\n");
84 
  la->A = PETSC_NULL;
85 
  la->B = PETSC_NULL;
86 
  ierr = STComputeExplicitOperator(eps->OP,&la->OP);CHKERRQ(ierr);
87 
  ierr = PetscTypeCompare((PetscObject)la->OP,MATSEQDENSE,&flg);CHKERRQ(ierr);
88 
  if (!flg) {
89 
    ierr = SlepcMatConvertSeqDense(la->OP,&la->OP);CHKERRQ(ierr);
90 
  }
1560 slepc 91 
  if (eps->extraction) {
92 
     ierr = PetscInfo(eps,"Warning: extraction type ignored\n");CHKERRQ(ierr);
1426 slepc 93 
  }
1596 slepc 94 
  ierr = EPSAllocateSolution(eps);CHKERRQ(ierr);
6 dsic.upv.es!jroman 95 
  PetscFunctionReturn(0);
96 
}
97 
 
98 
#undef __FUNCT__  
99 
#define __FUNCT__ "EPSSolve_LAPACK"
476 dsic.upv.es!antodo 100 
PetscErrorCode EPSSolve_LAPACK(EPS eps)
6 dsic.upv.es!jroman 101 
{
476 dsic.upv.es!antodo 102 
  PetscErrorCode ierr;
1928 jroman 103 
  PetscInt       n=eps->n,i,low,high;
982 slepc 104 
  PetscMPIInt    size;
780 dsic.upv.es!jroman 105 
  PetscScalar    *array,*arrayb,*pV,*pW;
679 dsic.upv.es!antodo 106 
  PetscReal      *w;
476 dsic.upv.es!antodo 107 
  EPS_LAPACK     *la = (EPS_LAPACK *)eps->data;
1422 slepc 108 
  MPI_Comm       comm = ((PetscObject)eps)->comm;
6 dsic.upv.es!jroman 109 
 
110 
  PetscFunctionBegin;
111 
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
393 dsic.upv.es!antodo 112 
  if (size == 1) {
113 
    ierr = VecGetArray(eps->V[0],&pV);CHKERRQ(ierr);
114 
  } else {
679 dsic.upv.es!antodo 115 
    ierr = PetscMalloc(sizeof(PetscScalar)*n*n,&pV);CHKERRQ(ierr);
393 dsic.upv.es!antodo 116 
  }
780 dsic.upv.es!jroman 117 
  if (eps->solverclass == EPS_TWO_SIDE && (la->OP || !eps->ishermitian)) {
118 
    if (size == 1) {
119 
      ierr = VecGetArray(eps->W[0],&pW);CHKERRQ(ierr);
120 
    } else {
121 
      ierr = PetscMalloc(sizeof(PetscScalar)*n*n,&pW);CHKERRQ(ierr);
122 
    }
123 
  } else pW = PETSC_NULL;
73 dsic.upv.es!antodo 124 
 
679 dsic.upv.es!antodo 125 
  if (la->OP) {
126 
    ierr = MatGetArray(la->OP,&array);CHKERRQ(ierr);
780 dsic.upv.es!jroman 127 
    ierr = EPSDenseNHEP(n,array,eps->eigr,eps->eigi,pV,pW);CHKERRQ(ierr);  
679 dsic.upv.es!antodo 128 
    ierr = MatRestoreArray(la->OP,&array);CHKERRQ(ierr);
129 
  } else if (eps->ishermitian) {
130 
#if defined(PETSC_USE_COMPLEX)
131 
    ierr = PetscMalloc(n*sizeof(PetscReal),&w);CHKERRQ(ierr);
132 
#else
133 
    w = eps->eigr;
134 
#endif
135 
    ierr = MatGetArray(la->A,&array);CHKERRQ(ierr);
136 
    if (!eps->isgeneralized) {
1175 slepc 137 
      ierr = EPSDenseHEP(n,array,n,w,pV);CHKERRQ(ierr);
679 dsic.upv.es!antodo 138 
    } else {
139 
      ierr = MatGetArray(la->B,&arrayb);CHKERRQ(ierr);
140 
      ierr = EPSDenseGHEP(n,array,arrayb,w,pV);CHKERRQ(ierr);  
141 
      ierr = MatRestoreArray(la->B,&arrayb);CHKERRQ(ierr);
142 
    }
143 
    ierr = MatRestoreArray(la->A,&array);CHKERRQ(ierr);
144 
#if defined(PETSC_USE_COMPLEX)
145 
    for (i=0;i<n;i++) {
146 
      eps->eigr[i] = w[i];
147 
    }
148 
    ierr = PetscFree(w);CHKERRQ(ierr);
149 
#endif    
150 
  } else {
151 
    ierr = MatGetArray(la->A,&array);CHKERRQ(ierr);
152 
    if (!eps->isgeneralized) {
780 dsic.upv.es!jroman 153 
      ierr = EPSDenseNHEP(n,array,eps->eigr,eps->eigi,pV,pW);CHKERRQ(ierr);
679 dsic.upv.es!antodo 154 
    } else {
155 
      ierr = MatGetArray(la->B,&arrayb);CHKERRQ(ierr);
780 dsic.upv.es!jroman 156 
      ierr = EPSDenseGNHEP(n,array,arrayb,eps->eigr,eps->eigi,pV,pW);CHKERRQ(ierr);  
679 dsic.upv.es!antodo 157 
      ierr = MatRestoreArray(la->B,&arrayb);CHKERRQ(ierr);
158 
    }
159 
    ierr = MatRestoreArray(la->A,&array);CHKERRQ(ierr);
160 
  }
6 dsic.upv.es!jroman 161 
 
393 dsic.upv.es!antodo 162 
  if (size == 1) {
163 
    ierr = VecRestoreArray(eps->V[0],&pV);CHKERRQ(ierr);
164 
  } else {
165 
    for (i=0; i<eps->ncv; i++) {
166 
      ierr = VecGetOwnershipRange(eps->V[i], &low, &high);CHKERRQ(ierr);
167 
      ierr = VecGetArray(eps->V[i], &array);CHKERRQ(ierr);
168 
      ierr = PetscMemcpy(array, pV+i*n+low, (high-low)*sizeof(PetscScalar));
169 
      ierr = VecRestoreArray(eps->V[i], &array);CHKERRQ(ierr);
73 dsic.upv.es!antodo 170 
    }
393 dsic.upv.es!antodo 171 
    ierr = PetscFree(pV);CHKERRQ(ierr);
6 dsic.upv.es!jroman 172 
  }
780 dsic.upv.es!jroman 173 
  if (pW) {
174 
    if (size == 1) {
175 
      ierr = VecRestoreArray(eps->W[0],&pW);CHKERRQ(ierr);
176 
    } else {
177 
      for (i=0; i<eps->ncv; i++) {
178 
        ierr = VecGetOwnershipRange(eps->W[i], &low, &high);CHKERRQ(ierr);
179 
        ierr = VecGetArray(eps->W[i], &array);CHKERRQ(ierr);
180 
        ierr = PetscMemcpy(array, pW+i*n+low, (high-low)*sizeof(PetscScalar));
181 
        ierr = VecRestoreArray(eps->W[i], &array);CHKERRQ(ierr);
182 
      }
183 
      ierr = PetscFree(pW);CHKERRQ(ierr);
184 
    }
185 
  }
6 dsic.upv.es!jroman 186 
 
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 
 
191 
  PetscFunctionReturn(0);
192 
}
193 
 
194 
#undef __FUNCT__  
195 
#define __FUNCT__ "EPSDestroy_LAPACK"
476 dsic.upv.es!antodo 196 
PetscErrorCode EPSDestroy_LAPACK(EPS eps)
6 dsic.upv.es!jroman 197 
{
476 dsic.upv.es!antodo 198 
  PetscErrorCode ierr;
199 
  EPS_LAPACK     *la = (EPS_LAPACK *)eps->data;
6 dsic.upv.es!jroman 200 
 
201 
  PetscFunctionBegin;
25 dsic.upv.es!jroman 202 
  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
679 dsic.upv.es!antodo 203 
  if (la->OP) { ierr = MatDestroy(la->OP);CHKERRQ(ierr); }
204 
  if (la->A) { ierr = MatDestroy(la->A);CHKERRQ(ierr); }
205 
  if (la->B) { ierr = MatDestroy(la->B);CHKERRQ(ierr); }
1040 slepc 206 
  ierr = PetscFree(eps->data);CHKERRQ(ierr);
1596 slepc 207 
  ierr = EPSFreeSolution(eps);CHKERRQ(ierr);
6 dsic.upv.es!jroman 208 
  PetscFunctionReturn(0);
209 
}
210 
 
211 
EXTERN_C_BEGIN
212 
#undef __FUNCT__  
213 
#define __FUNCT__ "EPSCreate_LAPACK"
476 dsic.upv.es!antodo 214 
PetscErrorCode EPSCreate_LAPACK(EPS eps)
6 dsic.upv.es!jroman 215 
{
476 dsic.upv.es!antodo 216 
  PetscErrorCode ierr;
217 
  EPS_LAPACK     *la;
6 dsic.upv.es!jroman 218 
 
219 
  PetscFunctionBegin;
220 
  ierr = PetscNew(EPS_LAPACK,&la);CHKERRQ(ierr);
221 
  PetscLogObjectMemory(eps,sizeof(EPS_LAPACK));
222 
  eps->data                      = (void *) la;
223 
  eps->ops->solve                = EPSSolve_LAPACK;
780 dsic.upv.es!jroman 224 
  eps->ops->solvets              = EPSSolve_LAPACK;
503 dsic.upv.es!antodo 225 
  eps->ops->setup                = EPSSetUp_LAPACK;
6 dsic.upv.es!jroman 226 
  eps->ops->destroy              = EPSDestroy_LAPACK;
177 dsic.upv.es!antodo 227 
  eps->ops->backtransform        = EPSBackTransform_Default;
503 dsic.upv.es!antodo 228 
  eps->ops->computevectors       = EPSComputeVectors_Default;
6 dsic.upv.es!jroman 229 
  PetscFunctionReturn(0);
230 
}
231 
EXTERN_C_END