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