Actual source code: samgpetsc.c
1: #define PETSCKSP_DLL
2: #include petscksp.h
3: #include src/mat/impls/aij/seq/aij.h
4: #include global.h
5: #include externc.h
6: #include samgfunc.h
7: #include petscfunc.h
9: /*MC
10: PCSAMG - SAMG + PETSc interface
11:
12: This interface allows e.g. to call samg as a shell preconditioner.
13: samg is the new algebraic multigrid code by Klaus Stueben [1].
14: Reference for SAMG
15: [1] K. St\"{u}ben,"Algebraic Multigrid: An Introduction for Positive
16: Definite Problems with Applications", in [2], pp. 413--532,
17: Academic Press, 2001.
18: [2] U. Trottenberg, C. Oosterlee and A. Sch\"{u}ller, "Multigrid
19: Methods", Academic Press, 2001.
20: [1] is also available as
21: [3] K. St\"{u}ben,"Algebraic Multigrid: An Introduction for Positive
22: Definite Problems with Applications", Tech. Rep. 53, German
23: National Research Center for Information Technology (GMD),
24: Schloss Birlinhoven, D-53754 Sankt-Augustin, Germany, March 1999
25: For more information on the SAMG-PETSc interface and examples of it's
26: use, see
27: [4] D. Lahaye "Algebraic Multigrid for Time-Harmonic Magnetic Field
28: Computations", PhD Thesis, KU Leuven, Belgium, December 2001.
29: (www.cs.kuleuven.ac.be/~domenico)
30:
31: Notes on PETSc part of this interface
32: Information on how to set up shell preconditioners in PETSc can be
33: in the PETSc documentation under preconditioners.
35: This preconditioner has not been completely organized to match the PETSc style,
36: see src/ksp/pc/impls/samgpetsc.c for the PC shell routines.
38: SAMG is a commercial product available from Klaus Stueben.
39:
40: Level: developer
42: Contributed by Domenico Lahaye (domenico.lahaye@cs.kuleuven.ac.be) January 2001
43:
44: M*/
46: /* ------------------------------------------------------------------- */
49: /*.. SamgShellPCCreate - This routine creates a user-defined
50: preconditioner context.
52: Output Parameter:
53: shell - user-defined preconditioner context..*/
55: PetscErrorCode SamgShellPCCreate(SamgShellPC **shell)
56: {
57: SamgShellPC *newctx;
60: PetscNew(SamgShellPC,&newctx);
61: *shell = newctx;
62: return 0;
63: }
65: /* ------------------------------------------------------------------- */
68: /*..SamgShellPCSetUp - This routine sets up a user-defined
69: ramg preconditioner context.
71: Input Parameters:
72: shell - user-defined preconditioner context
73: pmat - preconditioner matrix
75: Output Parameter:
76: shell - fully set up user-defined preconditioner context
78: This routine calls the setup phase of RAMG..*/
79:
80: PetscErrorCode SamgShellPCSetUp(SamgShellPC *shell, Mat pmat)
81: {
83: int numnodes, numnonzero, nnz_count;
84: int j, I, J, ncols_getrow, *cols_getrow, *diag;
85: PetscScalar *vals_getrow;
86: MatInfo info;
87: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)pmat->data;
88: /*..SAMG variables..*/
89: SAMG_PARAM *samg_param;
90: double *u_approx, *rhs, *Asky;
91: int *ia, *ja;
92: /*..Primary SAMG parameters..*/
93: /*....System of equations....*/
94: int matrix;
95: /*....Start solution and stopping criterium....*/
96: int ifirst;
97: double eps;
98: /*....Scalar and coupled system..*/
99: int nsys=1, ndiu=1, ndip=1;
100: int iscale[1], iu[1], ip[1];
101: /*....Approach and smoother....*/
102: int nsolve;
103: /*....Cycling process....*/
104: int ncyc;
105: /*....Repeated calls....*/
106: int iswtch;
107: /*....Initial dimensioning..*/
108: double a_cmplx, g_cmplx, p_cmplx, w_avrge;
109: /*....Class 1 input parameters controlling input/output..*/
110: int idump, iout;
111: double chktol;
112: /*....Output parameters....*/
113: double res_in, res_out;
114: int ncyc_done;
115: /*..Numbers of levels created by SAMG..*/
116: int levels;
117: /*..Auxilary integer to set secondary parameters..*/
118: int intin;
120: /*..Get size and number of unknowns of preconditioner matrix..*/
121: MatGetSize(pmat, &numnodes, &numnodes);
122: MatGetInfo(pmat,MAT_LOCAL,&info);
123: numnonzero = int(info.nz_used);
125: /*..Allocate memory for RAMG variables..*/
126: PetscMalloc(numnonzero * sizeof(double),&Asky);
127: PetscMalloc((numnodes+1) * sizeof(int),&ia);
128: PetscMalloc(numnonzero * sizeof(int),&ja);
129: PetscMalloc(numnodes * sizeof(double),&u_approx);
130: PetscMalloc(numnodes * sizeof(double),&rhs);
132: /*..Store PETSc matrix in compressed skyline format required by SAMG..*/
133: nnz_count = 0;
134: MatMarkDiagonal_SeqAIJ(pmat);
135: diag = aij->diag;
137: for (I=0;I<numnodes;I++){
138: ia[I] = nnz_count;
140: /*....put in diagonal entry first....*/
141: ja[nnz_count] = I;
142: Asky[nnz_count] = aij->a[diag[I]];
143: nnz_count++;
145: /*....put in off-diagonals....*/
146: ncols_getrow = aij->i[I+1] - aij->i[I];
147: vals_getrow = aij->a + aij->i[I];
148: cols_getrow = aij->j + aij->i[I];
149: for (j=0;j<ncols_getrow;j++){
150: J = cols_getrow[j];
151: if (J != I) {
152: Asky[nnz_count] = vals_getrow[j];
153: ja[nnz_count] = J;
154: nnz_count++;
155: }
156: }
157: }
158: ia[numnodes] = nnz_count;
160: /*..Allocate memory for SAMG parameters..*/
161: PetscNew(SAMG_PARAM,&samg_param);
163: /*..Set SAMG parameters..*/
164: SamgGetParam(samg_param);
166: /*..Set Primary parameters..*/
167: matrix = 12;
168: ifirst = (*samg_param).IFIRST;
169: eps = (*samg_param).EPS;
170: nsolve = (*samg_param).NSOLVE;
171: ncyc = (*samg_param).NCYC;
172: iswtch = (*samg_param).ISWTCH;
173: a_cmplx = (*samg_param).A_CMPLX;
174: g_cmplx = (*samg_param).G_CMPLX;
175: p_cmplx = (*samg_param).P_CMPLX;
176: w_avrge = (*samg_param).W_AVRGE;
177: chktol = (*samg_param).CHKTOL;
178: idump = (*samg_param).IDUMP;
179: iout = (*samg_param).IOUT;
180: /*..Set secondary parameters..*/
181: SAMG_set_levelx(&(samg_param->LEVELX));
182: SAMG_set_nptmn(&(samg_param->NPTMN));
183: SAMG_set_ecg(&(samg_param->ECG));
184: SAMG_set_ewt(&(samg_param->EWT));
185: SAMG_set_ncg(&(samg_param->NCG));
186: SAMG_set_nwt(&(samg_param->NWT));
187: SAMG_set_etr(&(samg_param->ETR));
188: SAMG_set_ntr(&(samg_param->NTR));
189: SAMG_set_nrd(&(samg_param->NRD));
190: SAMG_set_nru(&(samg_param->NRU));
191: SAMG_set_nrc(&(samg_param->NRC));
192: intin = 6; SAMG_set_logio(&intin);
193: intin = 1; SAMG_set_mode_debug(&intin);
194: intin = 0; SAMG_set_iter_pre(&intin);
195: intin = 0; SAMG_set_np_opt(&intin);
196: intin = 0; SAMG_set_ncgrad_default(&intin);
198: /*..Reset ncyc such that only setup is performed. This is done by setting
199: the last digit of ncyc (the number of cycles performed) equal to zero.
200: The first digits of ncyc (related to the solve phase) become
201: irrelevant.
202: ..*/
203: ncyc = 1000;
205: /*..Switch arrays ia and ja to Fortran conventions..*/
206: for (j=0;j<=numnodes;j++)
207: ia[j]++;
208: for (j=0;j<numnonzero;j++)
209: ja[j]++;
211: PetscPrintf(MPI_COMM_WORLD,"\n\n");
212: PetscPrintf(MPI_COMM_WORLD,"******************************************\n");
213: PetscPrintf(MPI_COMM_WORLD,"*** Start Setup SAMG code (scal. mode) ***\n");
214: PetscPrintf(MPI_COMM_WORLD,"******************************************\n");
215: PetscPrintf(MPI_COMM_WORLD,"\n\n");
217: /*..Call SAMG..*/
218: SAMG(&numnodes, &numnonzero, &nsys,
219: ia, ja, Asky, rhs, u_approx, iu, &ndiu, ip, &ndip, &matrix,
220: iscale, &res_in, &res_out, &ncyc_done, &ierr, &nsolve,
221: &ifirst, &eps, &ncyc, &iswtch, &a_cmplx, &g_cmplx,
222: &p_cmplx, &w_avrge, &chktol, &idump, &iout);
224: PetscPrintf(MPI_COMM_WORLD,"\n\n");
225: PetscPrintf(MPI_COMM_WORLD,"******************************************\n");
226: PetscPrintf(MPI_COMM_WORLD,"*** End Setup SAMG code (scal. mode) ***\n");
227: PetscPrintf(MPI_COMM_WORLD,"******************************************\n");
228: PetscPrintf(MPI_COMM_WORLD,"\n\n");
230: /*..Get number of levels created..*/
231: SAMGPETSC_get_levels(&levels);
233: /*..Store RAMG output in PETSc context..*/
234: shell->A = Asky;
235: shell->IA = ia;
236: shell->JA = ja;
237: shell->PARAM = samg_param;
238: (*shell).LEVELS = levels;
240: return 0;
241: }
243: /* ------------------------------------------------------------------- */
246: /*..SamgShellPCApply - This routine applies the AMG code as preconditioner
248: Input Parameters:
249: ctx - user-defined context, as set by SamgShellSetApply()
250: x - input vector
252: Output Parameter:
253: y - preconditioned vector
255: Notes:
256: Note that the PCSHELL preconditioner passes a void pointer as the
257: first input argument. This can be cast to be the whatever the user
258: has set (via PCSetShellApply()) the application-defined context to be.
260: ..*/
261: /*..To apply AMG as a preconditioner we set: */
262: /* i) rhs-vector equal to the residual */
263: /* ii) start solution equal to zero */
264: /* Implementation notes: */
265: /* For the residual (vector r) we take the values from the vector */
266: /* using VecGetArray. No explicit memory allocation for */
267: /* vals_getarray is thus needed as VecGetArray takes care of it. */
268: /* The allocated memory is freed again using a call to */
269: /* VecRestoreArray. The values in vals_getarray are then copy to */
270: /* rhs of the AMG code using memcpy. */
272: PetscErrorCode SamgShellPCApply(void *ctx, Vec r, Vec z)
273: {
275: int I, numnodes, numnonzero, *cols;
276: SamgShellPC *shell = (SamgShellPC *) ctx;
277: double *u_approx, *rhs, *Asky, *vals_getarray;
278: int *ia, *ja;
279: SAMG_PARAM *samg_param;
280: /*..Primary SAMG parameters..*/
281: /*....System of equations....*/
282: int matrix;
283: /*....Start solution and stopping criterium....*/
284: int ifirst;
285: double eps;
286: /*....Scalar and coupled system..*/
287: int nsys=1, ndiu=1, ndip=1;
288: int iscale[1], iu[1], ip[1];
289: /*....Approach and smoother....*/
290: int nsolve;
291: /*....Cycling process....*/
292: int ncyc;
293: /*....Repeated calls....*/
294: int iswtch;
295: /*....Initial dimensioning..*/
296: double a_cmplx, g_cmplx, p_cmplx, w_avrge;
297: /*....Class 1 input parameters controlling input/output..*/
298: int idump, iout;
299: double chktol;
300: /*....Output parameters....*/
301: double res_in, res_out;
302: int ncyc_done;
304: /*..Get values from context..*/
305: Asky = shell->A;
306: ia = shell->IA;
307: ja = shell->JA;
308: samg_param = shell->PARAM;
310: /*..Get numnodes and numnonzeros..*/
311: /*....numnodes can be determined as the size of the input vector r....*/
312: VecGetSize(r,&numnodes);
313: /*....numnonzero is determined from the pointer ia....*/
314: /*....Remember that ia following Fortran conventions....*/
315: numnonzero = ia[numnodes]-1;
317: /*..Set the rhs of the call to ramg equal to the residual..*/
318: VecGetArray(r,&vals_getarray);
320: /*..Allocate memory for rhs and initial solution of call to samg..*/
321: PetscMalloc(numnodes * sizeof(double),&u_approx);
322: PetscMalloc(numnodes * sizeof(double),&rhs);
324: /*..Set rhs of call to ramg..*/
325: memcpy(rhs, vals_getarray, numnodes * sizeof(*rhs));
326:
327: /*..Set initial solution of call to ramg to zero..*/
328: for (I=0;I<numnodes;I++){
329: u_approx[I] = 0.;
330: }
332: /*..Set Primary parameters..*/
333: matrix = (*samg_param).MATRIX;
334: ifirst = (*samg_param).IFIRST;
335: eps = (*samg_param).EPS;
336: nsolve = (*samg_param).NSOLVE;
337: ncyc = (*samg_param).NCYC;
338: iswtch = (*samg_param).ISWTCH;
339: a_cmplx = (*samg_param).A_CMPLX;
340: g_cmplx = (*samg_param).G_CMPLX;
341: p_cmplx = (*samg_param).P_CMPLX;
342: w_avrge = (*samg_param).W_AVRGE;
343: chktol = (*samg_param).CHKTOL;
344: idump = (*samg_param).IDUMP;
345: iout = (*samg_param).IOUT;
347: /*..Redefine iswtch to bypass setup..*/
348: /*....First digit of iswtch = 2: bypass setup and do not release memory
349: upon return....*/
350: /*....Second digit of iswtch = 1: memory extension switch....*/
351: /*....Third and fourth digit of iswtch: n_default. If n_default = 0,
352: the user has to set secondary parameters....*/
353: iswtch = 210;
355: /*..Call SAMG..*/
356: SAMG(&numnodes, &numnonzero, &nsys,
357: ia, ja, Asky, rhs, u_approx, iu, &ndiu, ip, &ndip, &matrix,
358: iscale, &res_in, &res_out, &ncyc_done, &ierr, &nsolve,
359: &ifirst, &eps, &ncyc, &iswtch, &a_cmplx, &g_cmplx,
360: &p_cmplx, &w_avrge, &chktol, &idump, &iout);
362: /*..Create auxilary integer array..*/
363: PetscMalloc(numnodes * sizeof(double),&cols);
365: for (I=0;I<numnodes;I++)
366: cols[I] = I;
368: /*..Store values computed by SAMG into the PETSc vector z..*/
369: VecSetValues(z,numnodes,cols,u_approx,INSERT_VALUES);
370:
372: /*..Restore PETSc rhs vector..*/
373: VecRestoreArray(r, &vals_getarray);
375: PetscFree(cols);
376: PetscFree(rhs);
377: PetscFree(u_approx);
379: return 0;
380: }
382: /* ------------------------------------------------------------------- */
385: /*..RamgShellPCDestroy - This routine destroys a user-defined
386: preconditioner context.
388: Input Parameter:
389: shell - user-defined preconditioner context..*/
391: PetscErrorCode SamgShellPCDestroy(SamgShellPC *shell)
392: {
393: /*..Free memory allocated by samg..*/
394: SAMG_cleanup();
395: /*..Free PCShell context..*/
396: PetscFree(shell->A);
397: PetscFree(shell->IA);
398: PetscFree(shell->JA);
399: PetscFree(shell->PARAM);
400: PetscFree(shell);
402: return 0;
403: }
404: /* ------------------------------------------------------------------- */
407: /*..SamgGetParam - Gets SAMG parameters specified at runtime
408: OUTPUT: The parameters set in the SAMG_PARAM context
409: ..*/
411: PetscErrorCode SamgGetParam(SAMG_PARAM *samg_param)
412: {
415: /*..Set default SAMG paramets..*/
416: /*....Class 0 SAMG parameters....*/
417: (*samg_param).MATRIX = 12;
418: /*....Primary SAMG parameters....*/
419: /*......If ifirst=0, the vector u, as passed to SAMG, is taken as first
420: approximation......*/
421: (*samg_param).IFIRST = 0;
422: (*samg_param).EPS = 1e-12;
423: /*......nsolve =1 denotes scalar approach......*/
424: (*samg_param).NSOLVE = 1;
425: /*......note: in the AMG-PETSc interface the number of cycles is required
426: to equal one to assure that in the PCApply routine AMG only performs
427: one cycle......*/
428: (*samg_param).NCYC = 1001;
429: (*samg_param).ISWTCH = 410;
430: (*samg_param).A_CMPLX = 2.5;
431: (*samg_param).G_CMPLX = 1.9;
432: (*samg_param).P_CMPLX = 1.9;
433: (*samg_param).W_AVRGE = 2.5;
434: (*samg_param).CHKTOL = 1e-8;
435: (*samg_param).IDUMP = 0;
436: (*samg_param).IOUT = 2;
437: /*....Secundary SAMG parameters....*/
438: (*samg_param).LEVELX = 25;
439: (*samg_param).NPTMN = 100;
440: (*samg_param).ECG = 0.25;
441: (*samg_param).EWT = 0.5;
442: (*samg_param).NCG = 1000;
443: (*samg_param).NWT = 3000;
444: (*samg_param).ETR = 12.2;
445: (*samg_param).NTR = 2;
446: (*samg_param).NRD = 131;
447: (*samg_param).NRC = 0;
448: (*samg_param).NRU = -131;
450: /*..Overwrite default values by values specified at runtime..*/
451: /*....Primary SAMG parameters....*/
452: PetscOptionsGetInt(PETSC_NULL,"-pc_samg_iswtch",&(*samg_param).ISWTCH,
453: PETSC_NULL);
455: PetscOptionsGetInt(PETSC_NULL,"-pc_samg_ncyc",&(*samg_param).NCYC,
456: PETSC_NULL);
458: PetscOptionsGetInt(PETSC_NULL,"-pc_samg_iout",&(*samg_param).IOUT,
459: PETSC_NULL);
461: /*....Secundary SAMG parameters....*/
462: PetscOptionsGetInt(PETSC_NULL,"-pc_samg_levelx",&(*samg_param).LEVELX,
463: PETSC_NULL);
465: PetscOptionsGetInt(PETSC_NULL,"-pc_samg_nptmn",&(*samg_param).NPTMN,
466: PETSC_NULL);
467: PetscOptionsGetInt(PETSC_NULL,"-pc_samg_nrd",&(*samg_param).NRD,
468: PETSC_NULL);
470: PetscOptionsGetInt(PETSC_NULL,"-pc_samg_nrc",&(*samg_param).NRC,
471: PETSC_NULL);
473: PetscOptionsGetInt(PETSC_NULL,"-pc_samg_nru",&(*samg_param).NRU,
474: PETSC_NULL);
476: return 0;
477: }