Actual source code: eige.c
1: /*$Id: eige.c,v 1.35 2001/08/07 03:03:45 balay Exp $*/
3: #include src/sles/ksp/kspimpl.h
5: #undef __FUNCT__
7: /*@
8: KSPComputeExplicitOperator - Computes the explicit preconditioned operator.
10: Collective on KSP
12: Input Parameter:
13: . ksp - the Krylov subspace context
15: Output Parameter:
16: . mat - the explict preconditioned operator
18: Notes:
19: This computation is done by applying the operators to columns of the
20: identity matrix.
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.
26: Level: advanced
27:
28: .keywords: KSP, compute, explicit, operator
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;
43: comm = ksp->comm;
45: MPI_Comm_size(comm,&size);
47: VecDuplicate(ksp->vec_sol,&in);
48: VecDuplicate(ksp->vec_sol,&out);
49: VecGetSize(in,&M);
50: VecGetLocalSize(in,&m);
51: VecGetOwnershipRange(in,&start,&end);
52: PetscMalloc((m+1)*sizeof(int),&rows);
53: for (i=0; i<m; i++) {rows[i] = start + i;}
55: if (size == 1) {
56: MatCreateSeqDense(comm,M,M,PETSC_NULL,mat);
57: } else {
58: /* MatCreateMPIDense(comm,m,M,M,M,PETSC_NULL,mat); */
59: MatCreateMPIAIJ(comm,m,m,M,M,0,0,0,0,mat);
60: }
61:
62: PCGetOperators(ksp->B,&A,PETSC_NULL,PETSC_NULL);
64: for (i=0; i<M; i++) {
66: VecSet(&zero,in);
67: VecSetValues(in,1,&i,&one,INSERT_VALUES);
68: VecAssemblyBegin(in);
69: VecAssemblyEnd(in);
71: KSP_MatMult(ksp,A,in,out);
72: KSP_PCApply(ksp,ksp->B,out,in);
73:
74: VecGetArray(in,&array);
75: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
76: VecRestoreArray(in,&array);
78: }
79: PetscFree(rows);
80: VecDestroy(in);
81: VecDestroy(out);
82: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
83: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
84: return(0);
85: }
87: #include petscblaslapack.h
89: #undef __FUNCT__
91: /*@
92: KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the
93: preconditioned operator using LAPACK.
95: Collective on KSP
97: Input Parameter:
98: + ksp - iterative context obtained from KSPCreate()
99: - n - size of arrays r and c
101: Output Parameters:
102: + r - real part of computed eigenvalues
103: - c - complex part of computed eigenvalues
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.
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.
115: Level: advanced
117: .keywords: KSP, compute, eigenvalues, explicitly
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;
132: KSPComputeExplicitOperator(ksp,&BA);
133: MPI_Comm_size(comm,&size);
134: MPI_Comm_rank(comm,&rank);
136: ierr = MatGetSize(BA,&n,&n);
137: if (size > 1) { /* assemble matrix on first processor */
138: if (!rank) {
139: MatCreateMPIDense(ksp->comm,n,n,n,n,PETSC_NULL,&A);
140: }
141: else {
142: MatCreateMPIDense(ksp->comm,0,n,n,n,PETSC_NULL,&A);
143: }
144: PetscLogObjectParent(BA,A);
146: MatGetOwnershipRange(BA,&row,&dummy);
147: MatGetLocalSize(BA,&m,&dummy);
148: for (i=0; i<m; i++) {
149: MatGetRow(BA,row,&nz,&cols,&vals);
150: MatSetValues(A,1,&row,nz,cols,vals,INSERT_VALUES);
151: MatRestoreRow(BA,row,&nz,&cols,&vals);
152: row++;
153: }
155: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
156: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
157: MatGetArray(A,&array);
158: } else {
159: MatGetArray(BA,&array);
160: }
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;
169: #if !defined(PETSC_USE_COMPLEX)
170: clen = n;
171: #else
172: clen = 2*n;
173: #endif
174: ierr = PetscMalloc(clen*sizeof(PetscScalar),&cwork);
175: idummy = n;
176: lwork = 5*n;
177: ierr = PetscMalloc(lwork*sizeof(PetscReal),&work);
178: ierr = PetscMalloc(n*sizeof(PetscReal),&realpart);
179: zero = 0;
180: LAgeev_(&zero,array,&n,cwork,&sdummy,&idummy,&idummy,&n,work,&lwork);
181: PetscFree(work);
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: PetscMalloc(n*sizeof(int),&perm);
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: PetscSortRealWithPermutation(n,realpart,perm);
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: PetscSortRealWithPermutation(n,realpart,perm);
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: PetscFree(perm);
209: PetscFree(realpart);
210: PetscFree(cwork);
211: }
212: #elif !defined(PETSC_USE_COMPLEX)
213: if (!rank) {
214: PetscScalar *work,sdummy;
215: PetscReal *realpart,*imagpart;
216: int idummy,lwork,*perm;
218: idummy = n;
219: lwork = 5*n;
220: ierr = PetscMalloc(2*n*sizeof(PetscReal),&realpart);
221: imagpart = realpart + n;
222: ierr = PetscMalloc(5*n*sizeof(PetscReal),&work);
223: #if defined(PETSC_MISSING_LAPACK_GEEV)
224: SETERRQ(PETSC_ERR_SUP,"GEEV - Lapack routine is unavilablenNot 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: PetscFree(work);
230: PetscMalloc(n*sizeof(int),&perm);
231: for (i=0; i<n; i++) { perm[i] = i;}
232: PetscSortRealWithPermutation(n,realpart,perm);
233: for (i=0; i<n; i++) {
234: r[i] = realpart[perm[i]];
235: c[i] = imagpart[perm[i]];
236: }
237: PetscFree(perm);
238: PetscFree(realpart);
239: }
240: #else
241: if (!rank) {
242: PetscScalar *work,sdummy,*eigs;
243: PetscReal *rwork;
244: int idummy,lwork,*perm;
246: idummy = n;
247: lwork = 5*n;
248: PetscMalloc(5*n*sizeof(PetscScalar),&work);
249: PetscMalloc(2*n*sizeof(PetscReal),&rwork);
250: PetscMalloc(n*sizeof(PetscScalar),&eigs);
251: #if defined(PETSC_MISSING_LAPACK_GEEV)
252: SETERRQ(PETSC_ERR_SUP,"GEEV - Lapack routine is unavilablenNot 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: PetscFree(work);
258: PetscFree(rwork);
259: PetscMalloc(n*sizeof(int),&perm);
260: for (i=0; i<n; i++) { perm[i] = i;}
261: for (i=0; i<n; i++) { r[i] = PetscRealPart(eigs[i]);}
262: PetscSortRealWithPermutation(n,r,perm);
263: for (i=0; i<n; i++) {
264: r[i] = PetscRealPart(eigs[perm[i]]);
265: c[i] = PetscImaginaryPart(eigs[perm[i]]);
266: }
267: PetscFree(perm);
268: PetscFree(eigs);
269: }
270: #endif
271: if (size > 1) {
272: MatRestoreArray(A,&array);
273: MatDestroy(A);
274: } else {
275: MatRestoreArray(BA,&array);
276: }
277: MatDestroy(BA);
278: return(0);
279: }