Actual source code: eige.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/kspimpl.h> /*I "petscksp.h" I*/
3: #include <petscblaslapack.h>
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(), PCComputeExplicitOperator()
31: @*/
32: PetscErrorCode KSPComputeExplicitOperator(KSP ksp,Mat *mat)
33: {
34: Vec in,out;
36: PetscMPIInt size;
37: PetscInt i,M,m,*rows,start,end;
38: Mat A;
39: MPI_Comm comm;
40: PetscScalar *array,one = 1.0;
45: comm = ((PetscObject)ksp)->comm;
47: MPI_Comm_size(comm,&size);
49: VecDuplicate(ksp->vec_sol,&in);
50: VecDuplicate(ksp->vec_sol,&out);
51: VecGetSize(in,&M);
52: VecGetLocalSize(in,&m);
53: VecGetOwnershipRange(in,&start,&end);
54: PetscMalloc(m*sizeof(PetscInt),&rows);
55: for (i=0; i<m; i++) {rows[i] = start + i;}
57: MatCreate(comm,mat);
58: MatSetSizes(*mat,m,m,M,M);
59: if (size == 1) {
60: MatSetType(*mat,MATSEQDENSE);
61: MatSeqDenseSetPreallocation(*mat,PETSC_NULL);
62: } else {
63: MatSetType(*mat,MATMPIAIJ);
64: MatMPIAIJSetPreallocation(*mat,0,PETSC_NULL,0,PETSC_NULL);
65: }
66: MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
67: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
68: PCGetOperators(ksp->pc,&A,PETSC_NULL,PETSC_NULL);
70: for (i=0; i<M; i++) {
72: VecSet(in,0.0);
73: VecSetValues(in,1,&i,&one,INSERT_VALUES);
74: VecAssemblyBegin(in);
75: VecAssemblyEnd(in);
77: KSP_MatMult(ksp,A,in,out);
78: KSP_PCApply(ksp,out,in);
79:
80: VecGetArray(in,&array);
81: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
82: VecRestoreArray(in,&array);
84: }
85: PetscFree(rows);
86: VecDestroy(&in);
87: VecDestroy(&out);
88: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
89: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
90: return(0);
91: }
95: /*@
96: KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the
97: preconditioned operator using LAPACK.
99: Collective on KSP
101: Input Parameter:
102: + ksp - iterative context obtained from KSPCreate()
103: - n - size of arrays r and c
105: Output Parameters:
106: + r - real part of computed eigenvalues
107: - c - complex part of computed eigenvalues
109: Notes:
110: This approach is very slow but will generally provide accurate eigenvalue
111: estimates. This routine explicitly forms a dense matrix representing
112: the preconditioned operator, and thus will run only for relatively small
113: problems, say n < 500.
115: Many users may just want to use the monitoring routine
116: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
117: to print the singular values at each iteration of the linear solve.
119: The preconditoner operator, rhs vector, solution vectors should be
120: set before this routine is called. i.e use KSPSetOperators(),KSPSolve() or
121: KSPSetOperators()
123: Level: advanced
125: .keywords: KSP, compute, eigenvalues, explicitly
127: .seealso: KSPComputeEigenvalues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues(), KSPSetOperators(), KSPSolve()
128: @*/
129: PetscErrorCode KSPComputeEigenvaluesExplicitly(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c)
130: {
131: Mat BA;
132: PetscErrorCode ierr;
133: PetscMPIInt size,rank;
134: MPI_Comm comm = ((PetscObject)ksp)->comm;
135: PetscScalar *array;
136: Mat A;
137: PetscInt m,row,nz,i,n,dummy;
138: const PetscInt *cols;
139: const PetscScalar *vals;
142: KSPComputeExplicitOperator(ksp,&BA);
143: MPI_Comm_size(comm,&size);
144: MPI_Comm_rank(comm,&rank);
146: MatGetSize(BA,&n,&n);
147: if (size > 1) { /* assemble matrix on first processor */
148: MatCreate(((PetscObject)ksp)->comm,&A);
149: if (!rank) {
150: MatSetSizes(A,n,n,n,n);
151: } else {
152: MatSetSizes(A,0,0,n,n);
153: }
154: MatSetType(A,MATMPIDENSE);
155: MatMPIDenseSetPreallocation(A,PETSC_NULL);
156: PetscLogObjectParent(BA,A);
158: MatGetOwnershipRange(BA,&row,&dummy);
159: MatGetLocalSize(BA,&m,&dummy);
160: for (i=0; i<m; i++) {
161: MatGetRow(BA,row,&nz,&cols,&vals);
162: MatSetValues(A,1,&row,nz,cols,vals,INSERT_VALUES);
163: MatRestoreRow(BA,row,&nz,&cols,&vals);
164: row++;
165: }
167: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
168: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
169: MatGetArray(A,&array);
170: } else {
171: MatGetArray(BA,&array);
172: }
174: #if defined(PETSC_HAVE_ESSL)
175: /* ESSL has a different calling sequence for dgeev() and zgeev() than standard LAPACK */
176: if (!rank) {
177: PetscScalar sdummy,*cwork;
178: PetscReal *work,*realpart;
179: PetscBLASInt clen,idummy,lwork,bn,zero = 0;
180: PetscInt *perm;
182: #if !defined(PETSC_USE_COMPLEX)
183: clen = n;
184: #else
185: clen = 2*n;
186: #endif
187: PetscMalloc(clen*sizeof(PetscScalar),&cwork);
188: idummy = -1; /* unused */
189: bn = PetscBLASIntCast(n);
190: lwork = 5*n;
191: PetscMalloc(lwork*sizeof(PetscReal),&work);
192: PetscMalloc(n*sizeof(PetscReal),&realpart);
193: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
194: LAPACKgeev_(&zero,array,&bn,cwork,&sdummy,&idummy,&idummy,&bn,work,&lwork);
195: PetscFPTrapPop();
196: PetscFree(work);
198: /* For now we stick with the convention of storing the real and imaginary
199: components of evalues separately. But is this what we really want? */
200: PetscMalloc(n*sizeof(PetscInt),&perm);
202: #if !defined(PETSC_USE_COMPLEX)
203: for (i=0; i<n; i++) {
204: realpart[i] = cwork[2*i];
205: perm[i] = i;
206: }
207: PetscSortRealWithPermutation(n,realpart,perm);
208: for (i=0; i<n; i++) {
209: r[i] = cwork[2*perm[i]];
210: c[i] = cwork[2*perm[i]+1];
211: }
212: #else
213: for (i=0; i<n; i++) {
214: realpart[i] = PetscRealPart(cwork[i]);
215: perm[i] = i;
216: }
217: PetscSortRealWithPermutation(n,realpart,perm);
218: for (i=0; i<n; i++) {
219: r[i] = PetscRealPart(cwork[perm[i]]);
220: c[i] = PetscImaginaryPart(cwork[perm[i]]);
221: }
222: #endif
223: PetscFree(perm);
224: PetscFree(realpart);
225: PetscFree(cwork);
226: }
227: #elif !defined(PETSC_USE_COMPLEX)
228: if (!rank) {
229: PetscScalar *work;
230: PetscReal *realpart,*imagpart;
231: PetscBLASInt idummy,lwork;
232: PetscInt *perm;
234: idummy = n;
235: lwork = 5*n;
236: PetscMalloc(2*n*sizeof(PetscReal),&realpart);
237: imagpart = realpart + n;
238: PetscMalloc(5*n*sizeof(PetscReal),&work);
239: #if defined(PETSC_MISSING_LAPACK_GEEV)
240: SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
241: #else
242: {
243: PetscBLASInt lierr;
244: PetscScalar sdummy;
245: PetscBLASInt bn = PetscBLASIntCast(n);
246: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
247: LAPACKgeev_("N","N",&bn,array,&bn,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr);
248: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
249: PetscFPTrapPop();
250: }
251: #endif
252: PetscFree(work);
253: PetscMalloc(n*sizeof(PetscInt),&perm);
254: for (i=0; i<n; i++) { perm[i] = i;}
255: PetscSortRealWithPermutation(n,realpart,perm);
256: for (i=0; i<n; i++) {
257: r[i] = realpart[perm[i]];
258: c[i] = imagpart[perm[i]];
259: }
260: PetscFree(perm);
261: PetscFree(realpart);
262: }
263: #else
264: if (!rank) {
265: PetscScalar *work,*eigs;
266: PetscReal *rwork;
267: PetscBLASInt idummy,lwork;
268: PetscInt *perm;
270: idummy = n;
271: lwork = 5*n;
272: PetscMalloc(5*n*sizeof(PetscScalar),&work);
273: PetscMalloc(2*n*sizeof(PetscReal),&rwork);
274: PetscMalloc(n*sizeof(PetscScalar),&eigs);
275: #if defined(PETSC_MISSING_LAPACK_GEEV)
276: SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
277: #else
278: {
279: PetscBLASInt lierr;
280: PetscScalar sdummy;
281: PetscBLASInt nb = PetscBLASIntCast(n);
282: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
283: LAPACKgeev_("N","N",&nb,array,&nb,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,rwork,&lierr);
284: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
285: PetscFPTrapPop();
286: }
287: #endif
288: PetscFree(work);
289: PetscFree(rwork);
290: PetscMalloc(n*sizeof(PetscInt),&perm);
291: for (i=0; i<n; i++) { perm[i] = i;}
292: for (i=0; i<n; i++) { r[i] = PetscRealPart(eigs[i]);}
293: PetscSortRealWithPermutation(n,r,perm);
294: for (i=0; i<n; i++) {
295: r[i] = PetscRealPart(eigs[perm[i]]);
296: c[i] = PetscImaginaryPart(eigs[perm[i]]);
297: }
298: PetscFree(perm);
299: PetscFree(eigs);
300: }
301: #endif
302: if (size > 1) {
303: MatRestoreArray(A,&array);
304: MatDestroy(&A);
305: } else {
306: MatRestoreArray(BA,&array);
307: }
308: MatDestroy(&BA);
309: return(0);
310: }