← Back to branch summary

~ubuntu-branches/ubuntu/feisty/petsc/feisty

« back to all changes in this revision

Viewing changes to src/sles/ksp/interface/eige.c

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2005-03-24 09:46:23 UTC
  • mfrom: (2.1.2 hoary)
  • Revision ID: james.westby@ubuntu.com-20050324094623-dfcxn8bltjms2cqq
Tags: 2.2.0-4
http://bugs.debian.org/284023
http://bugs.debian.org/293011
http://bugs.debian.org/274009
http://bugs.debian.org/267796
* Update for new mpich >> 1.2.5.3-2.
* Fixed src/inline/axpy.h for complex and UNROLL (closes: #284023, #293011).
* Added -fno-strict-aliasing to compile flags (closes: #274009).
* Switched SLES stuff to KSP in petsc.m4 (closes: 267796).

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/sles/ksp/interface/eige.c
1
 
/*$Id: eige.c,v 1.35 2001/08/07 03:03:45 balay Exp $*/
2
 
 
3
 
#include "src/sles/ksp/kspimpl.h"   /*I "petscksp.h" I*/
4
 
 
5
 
#undef __FUNCT__  
6
 
#define __FUNCT__ "KSPComputeExplicitOperator"
7
 
/*@
8
 
    KSPComputeExplicitOperator - Computes the explicit preconditioned operator.  
9
 
 
10
 
    Collective on KSP
11
 
 
12
 
    Input Parameter:
13
 
.   ksp - the Krylov subspace context
14
 
 
15
 
    Output Parameter:
16
 
.   mat - the explict preconditioned operator
17
 
 
18
 
    Notes:
19
 
    This computation is done by applying the operators to columns of the 
20
 
    identity matrix.
21
 
 
22
 
    Currently, this routine uses a dense matrix format when 1 processor
23
 
    is used and a sparse format otherwise.  This routine is costly in general,
24
 
    and is recommended for use only with relatively small systems.
25
 
 
26
 
    Level: advanced
27
 
   
28
 
.keywords: KSP, compute, explicit, operator
29
 
 
30
 
.seealso: KSPComputeEigenvaluesExplicitly()
31
 
@*/
32
 
int KSPComputeExplicitOperator(KSP ksp,Mat *mat)
33
 
{
34
 
  Vec      in,out;
35
 
  int      ierr,i,M,m,size,*rows,start,end;
36
 
  Mat      A;
37
 
  MPI_Comm comm;
38
 
  PetscScalar   *array,zero = 0.0,one = 1.0;
39
 
 
40
 
  PetscFunctionBegin;
41
 
  PetscValidHeaderSpecific(ksp,KSP_COOKIE);
42
 
  PetscValidPointer(mat);
43
 
  comm = ksp->comm;
44
 
 
45
 
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
46
 
 
47
 
  ierr = VecDuplicate(ksp->vec_sol,&in);CHKERRQ(ierr);
48
 
  ierr = VecDuplicate(ksp->vec_sol,&out);CHKERRQ(ierr);
49
 
  ierr = VecGetSize(in,&M);CHKERRQ(ierr);
50
 
  ierr = VecGetLocalSize(in,&m);CHKERRQ(ierr);
51
 
  ierr = VecGetOwnershipRange(in,&start,&end);CHKERRQ(ierr);
52
 
  ierr = PetscMalloc((m+1)*sizeof(int),&rows);CHKERRQ(ierr);
53
 
  for (i=0; i<m; i++) {rows[i] = start + i;}
54
 
 
55
 
  if (size == 1) {
56
 
    ierr = MatCreateSeqDense(comm,M,M,PETSC_NULL,mat);CHKERRQ(ierr);
57
 
  } else {
58
 
    /*    ierr = MatCreateMPIDense(comm,m,M,M,M,PETSC_NULL,mat);CHKERRQ(ierr); */
59
 
    ierr = MatCreateMPIAIJ(comm,m,m,M,M,0,0,0,0,mat);CHKERRQ(ierr);
60
 
  }
61
 
  
62
 
  ierr = PCGetOperators(ksp->B,&A,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
63
 
 
64
 
  for (i=0; i<M; i++) {
65
 
 
66
 
    ierr = VecSet(&zero,in);CHKERRQ(ierr);
67
 
    ierr = VecSetValues(in,1,&i,&one,INSERT_VALUES);CHKERRQ(ierr);
68
 
    ierr = VecAssemblyBegin(in);CHKERRQ(ierr);
69
 
    ierr = VecAssemblyEnd(in);CHKERRQ(ierr);
70
 
 
71
 
    ierr = KSP_MatMult(ksp,A,in,out);CHKERRQ(ierr);
72
 
    ierr = KSP_PCApply(ksp,ksp->B,out,in);CHKERRQ(ierr);
73
 
    
74
 
    ierr = VecGetArray(in,&array);CHKERRQ(ierr);
75
 
    ierr = MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 
76
 
    ierr = VecRestoreArray(in,&array);CHKERRQ(ierr);
77
 
 
78
 
  }
79
 
  ierr = PetscFree(rows);CHKERRQ(ierr);
80
 
  ierr = VecDestroy(in);CHKERRQ(ierr);
81
 
  ierr = VecDestroy(out);CHKERRQ(ierr);
82
 
  ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
83
 
  ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
84
 
  PetscFunctionReturn(0);
85
 
}
86
 
 
87
 
#include "petscblaslapack.h"
88
 
 
89
 
#undef __FUNCT__  
90
 
#define __FUNCT__ "KSPComputeEigenvaluesExplicitly"
91
 
/*@
92
 
   KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the 
93
 
   preconditioned operator using LAPACK.  
94
 
 
95
 
   Collective on KSP
96
 
 
97
 
   Input Parameter:
98
 
+  ksp - iterative context obtained from KSPCreate()
99
 
-  n - size of arrays r and c
100
 
 
101
 
   Output Parameters:
102
 
+  r - real part of computed eigenvalues
103
 
-  c - complex part of computed eigenvalues
104
 
 
105
 
   Notes:
106
 
   This approach is very slow but will generally provide accurate eigenvalue
107
 
   estimates.  This routine explicitly forms a dense matrix representing 
108
 
   the preconditioned operator, and thus will run only for relatively small
109
 
   problems, say n < 500.
110
 
 
111
 
   Many users may just want to use the monitoring routine
112
 
   KSPSingularValueMonitor() (which can be set with option -ksp_singmonitor)
113
 
   to print the singular values at each iteration of the linear solve.
114
 
 
115
 
   Level: advanced
116
 
 
117
 
.keywords: KSP, compute, eigenvalues, explicitly
118
 
 
119
 
.seealso: KSPComputeEigenvalues(), KSPSingularValueMonitor(), KSPComputeExtremeSingularValues()
120
 
@*/
121
 
int KSPComputeEigenvaluesExplicitly(KSP ksp,int nmax,PetscReal *r,PetscReal *c) 
122
 
{
123
 
  Mat          BA;
124
 
  int          i,n,ierr,size,rank,dummy;
125
 
  MPI_Comm     comm = ksp->comm;
126
 
  PetscScalar  *array;
127
 
  Mat          A;
128
 
  int          m,row,nz,*cols;
129
 
  PetscScalar  *vals;
130
 
 
131
 
  PetscFunctionBegin;
132
 
  ierr =  KSPComputeExplicitOperator(ksp,&BA);CHKERRQ(ierr);
133
 
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
134
 
  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
135
 
 
136
 
  ierr     = MatGetSize(BA,&n,&n);CHKERRQ(ierr);
137
 
  if (size > 1) { /* assemble matrix on first processor */
138
 
    if (!rank) {
139
 
      ierr = MatCreateMPIDense(ksp->comm,n,n,n,n,PETSC_NULL,&A);CHKERRQ(ierr);
140
 
    }
141
 
    else {
142
 
      ierr = MatCreateMPIDense(ksp->comm,0,n,n,n,PETSC_NULL,&A);CHKERRQ(ierr);
143
 
    }
144
 
    PetscLogObjectParent(BA,A);
145
 
 
146
 
    ierr = MatGetOwnershipRange(BA,&row,&dummy);CHKERRQ(ierr);
147
 
    ierr = MatGetLocalSize(BA,&m,&dummy);CHKERRQ(ierr);
148
 
    for (i=0; i<m; i++) {
149
 
      ierr = MatGetRow(BA,row,&nz,&cols,&vals);CHKERRQ(ierr);
150
 
      ierr = MatSetValues(A,1,&row,nz,cols,vals,INSERT_VALUES);CHKERRQ(ierr);
151
 
      ierr = MatRestoreRow(BA,row,&nz,&cols,&vals);CHKERRQ(ierr);
152
 
      row++;
153
 
    } 
154
 
 
155
 
    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
156
 
    ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
157
 
    ierr = MatGetArray(A,&array);CHKERRQ(ierr);
158
 
  } else {
159
 
    ierr = MatGetArray(BA,&array);CHKERRQ(ierr);
160
 
  }
161
 
 
162
 
#if defined(PETSC_HAVE_ESSL)
163
 
  /* ESSL has a different calling sequence for dgeev() and zgeev() than standard LAPACK */
164
 
  if (!rank) {
165
 
    PetscScalar    sdummy,*cwork;
166
 
    PetscReal *work,*realpart;
167
 
    int       clen,idummy,lwork,*perm,zero;
168
 
 
169
 
#if !defined(PETSC_USE_COMPLEX)
170
 
    clen = n;
171
 
#else
172
 
    clen = 2*n;
173
 
#endif
174
 
    ierr   = PetscMalloc(clen*sizeof(PetscScalar),&cwork);CHKERRQ(ierr);
175
 
    idummy = n;
176
 
    lwork  = 5*n;
177
 
    ierr   = PetscMalloc(lwork*sizeof(PetscReal),&work);CHKERRQ(ierr);
178
 
    ierr   = PetscMalloc(n*sizeof(PetscReal),&realpart);CHKERRQ(ierr);
179
 
    zero   = 0;
180
 
    LAgeev_(&zero,array,&n,cwork,&sdummy,&idummy,&idummy,&n,work,&lwork);
181
 
    ierr = PetscFree(work);CHKERRQ(ierr);
182
 
 
183
 
    /* For now we stick with the convention of storing the real and imaginary
184
 
       components of evalues separately.  But is this what we really want? */
185
 
    ierr = PetscMalloc(n*sizeof(int),&perm);CHKERRQ(ierr);
186
 
 
187
 
#if !defined(PETSC_USE_COMPLEX)
188
 
    for (i=0; i<n; i++) {
189
 
      realpart[i] = cwork[2*i];
190
 
      perm[i]     = i;
191
 
    }
192
 
    ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr);
193
 
    for (i=0; i<n; i++) {
194
 
      r[i] = cwork[2*perm[i]];
195
 
      c[i] = cwork[2*perm[i]+1];
196
 
    }
197
 
#else
198
 
    for (i=0; i<n; i++) {
199
 
      realpart[i] = PetscRealPart(cwork[i]);
200
 
      perm[i]     = i;
201
 
    }
202
 
    ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr);
203
 
    for (i=0; i<n; i++) {
204
 
      r[i] = PetscRealPart(cwork[perm[i]]);
205
 
      c[i] = PetscImaginaryPart(cwork[perm[i]]);
206
 
    }
207
 
#endif
208
 
    ierr = PetscFree(perm);CHKERRQ(ierr);
209
 
    ierr = PetscFree(realpart);CHKERRQ(ierr);
210
 
    ierr = PetscFree(cwork);CHKERRQ(ierr);
211
 
  }
212
 
#elif !defined(PETSC_USE_COMPLEX)
213
 
  if (!rank) {
214
 
    PetscScalar    *work,sdummy;
215
 
    PetscReal *realpart,*imagpart;
216
 
    int       idummy,lwork,*perm;
217
 
 
218
 
    idummy   = n;
219
 
    lwork    = 5*n;
220
 
    ierr     = PetscMalloc(2*n*sizeof(PetscReal),&realpart);CHKERRQ(ierr);
221
 
    imagpart = realpart + n;
222
 
    ierr     = PetscMalloc(5*n*sizeof(PetscReal),&work);CHKERRQ(ierr);
223
 
#if defined(PETSC_MISSING_LAPACK_GEEV) 
224
 
  SETERRQ(PETSC_ERR_SUP,"GEEV - Lapack routine is unavilable\nNot able to provide eigen values.");
225
 
#else
226
 
    LAgeev_("N","N",&n,array,&n,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&ierr);
227
 
#endif
228
 
    if (ierr) SETERRQ1(PETSC_ERR_LIB,"Error in LAPACK routine %d",ierr);
229
 
    ierr = PetscFree(work);CHKERRQ(ierr);
230
 
    ierr = PetscMalloc(n*sizeof(int),&perm);CHKERRQ(ierr);
231
 
    for (i=0; i<n; i++) { perm[i] = i;}
232
 
    ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr);
233
 
    for (i=0; i<n; i++) {
234
 
      r[i] = realpart[perm[i]];
235
 
      c[i] = imagpart[perm[i]];
236
 
    }
237
 
    ierr = PetscFree(perm);CHKERRQ(ierr);
238
 
    ierr = PetscFree(realpart);CHKERRQ(ierr);
239
 
  }
240
 
#else
241
 
  if (!rank) {
242
 
    PetscScalar *work,sdummy,*eigs;
243
 
    PetscReal *rwork;
244
 
    int    idummy,lwork,*perm;
245
 
 
246
 
    idummy   = n;
247
 
    lwork    = 5*n;
248
 
    ierr = PetscMalloc(5*n*sizeof(PetscScalar),&work);CHKERRQ(ierr);
249
 
    ierr = PetscMalloc(2*n*sizeof(PetscReal),&rwork);CHKERRQ(ierr);
250
 
    ierr = PetscMalloc(n*sizeof(PetscScalar),&eigs);CHKERRQ(ierr);
251
 
#if defined(PETSC_MISSING_LAPACK_GEEV) 
252
 
  SETERRQ(PETSC_ERR_SUP,"GEEV - Lapack routine is unavilable\nNot able to provide eigen values.");
253
 
#else
254
 
    LAgeev_("N","N",&n,array,&n,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,rwork,&ierr);
255
 
#endif
256
 
    if (ierr) SETERRQ1(PETSC_ERR_LIB,"Error in LAPACK routine %d",ierr);
257
 
    ierr = PetscFree(work);CHKERRQ(ierr);
258
 
    ierr = PetscFree(rwork);CHKERRQ(ierr);
259
 
    ierr = PetscMalloc(n*sizeof(int),&perm);CHKERRQ(ierr);
260
 
    for (i=0; i<n; i++) { perm[i] = i;}
261
 
    for (i=0; i<n; i++) { r[i]    = PetscRealPart(eigs[i]);}
262
 
    ierr = PetscSortRealWithPermutation(n,r,perm);CHKERRQ(ierr);
263
 
    for (i=0; i<n; i++) {
264
 
      r[i] = PetscRealPart(eigs[perm[i]]);
265
 
      c[i] = PetscImaginaryPart(eigs[perm[i]]);
266
 
    }
267
 
    ierr = PetscFree(perm);CHKERRQ(ierr);
268
 
    ierr = PetscFree(eigs);CHKERRQ(ierr);
269
 
  }
270
 
#endif  
271
 
  if (size > 1) {
272
 
    ierr = MatRestoreArray(A,&array);CHKERRQ(ierr);
273
 
    ierr = MatDestroy(A);CHKERRQ(ierr);
274
 
  } else {
275
 
    ierr = MatRestoreArray(BA,&array);CHKERRQ(ierr);
276
 
  }
277
 
  ierr = MatDestroy(BA);CHKERRQ(ierr);
278
 
  PetscFunctionReturn(0);
279
 
}