Subversion Repositories slepc-dev

Rev

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

Rev Author Line No. Line
6 dsic.upv.es!jroman 1 
 
2 
/*                      
3 
       This file implements a wrapper to the LAPACK eigenvalue subroutines.
679 dsic.upv.es!antodo 4 
       Generalized problems are transformed to standard ones only if necessary.
6 dsic.upv.es!jroman 5 
*/
679 dsic.upv.es!antodo 6 
#include "src/eps/epsimpl.h"
6 dsic.upv.es!jroman 7 
#include "slepcblaslapack.h"
8 
 
679 dsic.upv.es!antodo 9 
typedef struct {
10 
  Mat OP,A,B;
11 
} EPS_LAPACK;
12 
 
6 dsic.upv.es!jroman 13 
#undef __FUNCT__  
14 
#define __FUNCT__ "EPSSetUp_LAPACK"
476 dsic.upv.es!antodo 15 
PetscErrorCode EPSSetUp_LAPACK(EPS eps)
6 dsic.upv.es!jroman 16 
{
476 dsic.upv.es!antodo 17 
  PetscErrorCode ierr;
982 slepc 18 
  PetscInt       N;
476 dsic.upv.es!antodo 19 
  EPS_LAPACK     *la = (EPS_LAPACK *)eps->data;
20 
  PetscTruth     flg;
679 dsic.upv.es!antodo 21 
  Mat            A,B;
22 
  PetscScalar    shift;
23 
 
6 dsic.upv.es!jroman 24 
  PetscFunctionBegin;
260 dsic.upv.es!antodo 25 
  ierr = VecGetSize(eps->vec_initial,&N);CHKERRQ(ierr);
679 dsic.upv.es!antodo 26 
  if (eps->nev<1 || eps->nev>N) SETERRQ(1,"Wrong value of nev");
27 
  eps->ncv = N;
260 dsic.upv.es!antodo 28 
 
679 dsic.upv.es!antodo 29 
  if (la->OP) { ierr = MatDestroy(la->OP);CHKERRQ(ierr); }
30 
  if (la->A) { ierr = MatDestroy(la->A);CHKERRQ(ierr); }
31 
  if (la->B) { ierr = MatDestroy(la->B);CHKERRQ(ierr); }
32 
 
33 
  ierr = PetscTypeCompare((PetscObject)eps->OP,STSHIFT,&flg);CHKERRQ(ierr);
34 
  if (flg) {
35 
    la->OP = PETSC_NULL;
36 
    ierr = STGetOperators(eps->OP,&A,&B);CHKERRQ(ierr);
37 
    ierr = SlepcMatConvertSeqDense(A,&la->A);CHKERRQ(ierr);
38 
    if (eps->isgeneralized) {
39 
      ierr = SlepcMatConvertSeqDense(B,&la->B);CHKERRQ(ierr);
40 
    } else la->B = PETSC_NULL;
41 
    ierr = STGetShift(eps->OP,&shift);CHKERRQ(ierr);
42 
    if (shift != 0.0) {
828 dsic.upv.es!antodo 43 
      ierr = MatShift(la->A,shift);CHKERRQ(ierr);
679 dsic.upv.es!antodo 44 
    }
45 
  } else {
1038 slepc 46 
    PetscInfo(eps,"Using slow explicit operator\n");
679 dsic.upv.es!antodo 47 
    la->A = PETSC_NULL;
48 
    la->B = PETSC_NULL;
49 
    ierr = STComputeExplicitOperator(eps->OP,&la->OP);CHKERRQ(ierr);
50 
    ierr = PetscTypeCompare((PetscObject)la->OP,MATSEQDENSE,&flg);CHKERRQ(ierr);
51 
    if (!flg) {
52 
      ierr = SlepcMatConvertSeqDense(la->OP,&la->OP);CHKERRQ(ierr);
53 
    }
6 dsic.upv.es!jroman 54 
  }
55 
 
260 dsic.upv.es!antodo 56 
  ierr = EPSAllocateSolutionContiguous(eps);CHKERRQ(ierr);
6 dsic.upv.es!jroman 57 
  PetscFunctionReturn(0);
58 
}
59 
 
60 
#undef __FUNCT__  
61 
#define __FUNCT__ "EPSSolve_LAPACK"
476 dsic.upv.es!antodo 62 
PetscErrorCode EPSSolve_LAPACK(EPS eps)
6 dsic.upv.es!jroman 63 
{
476 dsic.upv.es!antodo 64 
  PetscErrorCode ierr;
982 slepc 65 
  PetscInt       n,i,low,high;
66 
  PetscMPIInt    size;
780 dsic.upv.es!jroman 67 
  PetscScalar    *array,*arrayb,*pV,*pW;
679 dsic.upv.es!antodo 68 
  PetscReal      *w;
476 dsic.upv.es!antodo 69 
  EPS_LAPACK     *la = (EPS_LAPACK *)eps->data;
70 
  MPI_Comm       comm = eps->comm;
6 dsic.upv.es!jroman 71 
 
72 
  PetscFunctionBegin;
73 
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
73 dsic.upv.es!antodo 74 
 
679 dsic.upv.es!antodo 75 
  ierr = VecGetSize(eps->vec_initial,&n);CHKERRQ(ierr);
6 dsic.upv.es!jroman 76 
 
393 dsic.upv.es!antodo 77 
  if (size == 1) {
78 
    ierr = VecGetArray(eps->V[0],&pV);CHKERRQ(ierr);
79 
  } else {
679 dsic.upv.es!antodo 80 
    ierr = PetscMalloc(sizeof(PetscScalar)*n*n,&pV);CHKERRQ(ierr);
393 dsic.upv.es!antodo 81 
  }
780 dsic.upv.es!jroman 82 
  if (eps->solverclass == EPS_TWO_SIDE && (la->OP || !eps->ishermitian)) {
83 
    if (size == 1) {
84 
      ierr = VecGetArray(eps->W[0],&pW);CHKERRQ(ierr);
85 
    } else {
86 
      ierr = PetscMalloc(sizeof(PetscScalar)*n*n,&pW);CHKERRQ(ierr);
87 
    }
88 
  } else pW = PETSC_NULL;
73 dsic.upv.es!antodo 89 
 
780 dsic.upv.es!jroman 90 
 
679 dsic.upv.es!antodo 91 
  if (la->OP) {
92 
    ierr = MatGetArray(la->OP,&array);CHKERRQ(ierr);
780 dsic.upv.es!jroman 93 
    ierr = EPSDenseNHEP(n,array,eps->eigr,eps->eigi,pV,pW);CHKERRQ(ierr);  
679 dsic.upv.es!antodo 94 
    ierr = MatRestoreArray(la->OP,&array);CHKERRQ(ierr);
95 
  } else if (eps->ishermitian) {
96 
#if defined(PETSC_USE_COMPLEX)
97 
    ierr = PetscMalloc(n*sizeof(PetscReal),&w);CHKERRQ(ierr);
98 
#else
99 
    w = eps->eigr;
100 
#endif
101 
    ierr = MatGetArray(la->A,&array);CHKERRQ(ierr);
102 
    if (!eps->isgeneralized) {
103 
      ierr = EPSDenseHEP(n,array,w,pV);CHKERRQ(ierr);
104 
    } else {
105 
      ierr = MatGetArray(la->B,&arrayb);CHKERRQ(ierr);
106 
      ierr = EPSDenseGHEP(n,array,arrayb,w,pV);CHKERRQ(ierr);  
107 
      ierr = MatRestoreArray(la->B,&arrayb);CHKERRQ(ierr);
108 
    }
109 
    ierr = MatRestoreArray(la->A,&array);CHKERRQ(ierr);
110 
#if defined(PETSC_USE_COMPLEX)
111 
    for (i=0;i<n;i++) {
112 
      eps->eigr[i] = w[i];
113 
    }
114 
    ierr = PetscFree(w);CHKERRQ(ierr);
115 
#endif    
683 dsic.upv.es!jroman 116 
    for (i=0;i<n;i++) eps->eigi[i]=0.0;
679 dsic.upv.es!antodo 117 
  } else {
118 
    ierr = MatGetArray(la->A,&array);CHKERRQ(ierr);
119 
    if (!eps->isgeneralized) {
780 dsic.upv.es!jroman 120 
      ierr = EPSDenseNHEP(n,array,eps->eigr,eps->eigi,pV,pW);CHKERRQ(ierr);
679 dsic.upv.es!antodo 121 
    } else {
122 
      ierr = MatGetArray(la->B,&arrayb);CHKERRQ(ierr);
780 dsic.upv.es!jroman 123 
      ierr = EPSDenseGNHEP(n,array,arrayb,eps->eigr,eps->eigi,pV,pW);CHKERRQ(ierr);  
679 dsic.upv.es!antodo 124 
      ierr = MatRestoreArray(la->B,&arrayb);CHKERRQ(ierr);
125 
    }
126 
    ierr = MatRestoreArray(la->A,&array);CHKERRQ(ierr);
127 
  }
6 dsic.upv.es!jroman 128 
 
393 dsic.upv.es!antodo 129 
  if (size == 1) {
130 
    ierr = VecRestoreArray(eps->V[0],&pV);CHKERRQ(ierr);
131 
  } else {
132 
    for (i=0; i<eps->ncv; i++) {
133 
      ierr = VecGetOwnershipRange(eps->V[i], &low, &high);CHKERRQ(ierr);
134 
      ierr = VecGetArray(eps->V[i], &array);CHKERRQ(ierr);
135 
      ierr = PetscMemcpy(array, pV+i*n+low, (high-low)*sizeof(PetscScalar));
136 
      ierr = VecRestoreArray(eps->V[i], &array);CHKERRQ(ierr);
73 dsic.upv.es!antodo 137 
    }
393 dsic.upv.es!antodo 138 
    ierr = PetscFree(pV);CHKERRQ(ierr);
6 dsic.upv.es!jroman 139 
  }
780 dsic.upv.es!jroman 140 
  if (pW) {
141 
    if (size == 1) {
142 
      ierr = VecRestoreArray(eps->W[0],&pW);CHKERRQ(ierr);
143 
    } else {
144 
      for (i=0; i<eps->ncv; i++) {
145 
        ierr = VecGetOwnershipRange(eps->W[i], &low, &high);CHKERRQ(ierr);
146 
        ierr = VecGetArray(eps->W[i], &array);CHKERRQ(ierr);
147 
        ierr = PetscMemcpy(array, pW+i*n+low, (high-low)*sizeof(PetscScalar));
148 
        ierr = VecRestoreArray(eps->W[i], &array);CHKERRQ(ierr);
149 
      }
150 
      ierr = PetscFree(pW);CHKERRQ(ierr);
151 
    }
152 
  }
6 dsic.upv.es!jroman 153 
 
154 
  eps->nconv = eps->ncv;
146 dsic.upv.es!antodo 155 
  eps->its   = 1;  
156 
  eps->reason = EPS_CONVERGED_TOL;
6 dsic.upv.es!jroman 157 
 
158 
  PetscFunctionReturn(0);
159 
}
160 
 
161 
#undef __FUNCT__  
162 
#define __FUNCT__ "EPSDestroy_LAPACK"
476 dsic.upv.es!antodo 163 
PetscErrorCode EPSDestroy_LAPACK(EPS eps)
6 dsic.upv.es!jroman 164 
{
476 dsic.upv.es!antodo 165 
  PetscErrorCode ierr;
166 
  EPS_LAPACK     *la = (EPS_LAPACK *)eps->data;
6 dsic.upv.es!jroman 167 
 
168 
  PetscFunctionBegin;
25 dsic.upv.es!jroman 169 
  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
679 dsic.upv.es!antodo 170 
  if (la->OP) { ierr = MatDestroy(la->OP);CHKERRQ(ierr); }
171 
  if (la->A) { ierr = MatDestroy(la->A);CHKERRQ(ierr); }
172 
  if (la->B) { ierr = MatDestroy(la->B);CHKERRQ(ierr); }
1040 slepc 173 
  ierr = PetscFree(eps->data);CHKERRQ(ierr);
260 dsic.upv.es!antodo 174 
  ierr = EPSFreeSolutionContiguous(eps);CHKERRQ(ierr);
6 dsic.upv.es!jroman 175 
  PetscFunctionReturn(0);
176 
}
177 
 
178 
EXTERN_C_BEGIN
179 
#undef __FUNCT__  
180 
#define __FUNCT__ "EPSCreate_LAPACK"
476 dsic.upv.es!antodo 181 
PetscErrorCode EPSCreate_LAPACK(EPS eps)
6 dsic.upv.es!jroman 182 
{
476 dsic.upv.es!antodo 183 
  PetscErrorCode ierr;
184 
  EPS_LAPACK     *la;
6 dsic.upv.es!jroman 185 
 
186 
  PetscFunctionBegin;
187 
  ierr = PetscNew(EPS_LAPACK,&la);CHKERRQ(ierr);
188 
  PetscLogObjectMemory(eps,sizeof(EPS_LAPACK));
189 
  eps->data                      = (void *) la;
190 
  eps->ops->solve                = EPSSolve_LAPACK;
780 dsic.upv.es!jroman 191 
  eps->ops->solvets              = EPSSolve_LAPACK;
503 dsic.upv.es!antodo 192 
  eps->ops->setup                = EPSSetUp_LAPACK;
6 dsic.upv.es!jroman 193 
  eps->ops->destroy              = EPSDestroy_LAPACK;
177 dsic.upv.es!antodo 194 
  eps->ops->backtransform        = EPSBackTransform_Default;
503 dsic.upv.es!antodo 195 
  eps->ops->computevectors       = EPSComputeVectors_Default;
6 dsic.upv.es!jroman 196 
  PetscFunctionReturn(0);
197 
}
198 
EXTERN_C_END