/[escript]/branches/escript3047_with_pastix2995/paso/src/PASTIX.c_debug
ViewVC logotype

Annotation of /branches/escript3047_with_pastix2995/paso/src/PASTIX.c_debug

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3049 - (hide annotations)
Fri Jun 25 04:20:29 2010 UTC (8 years, 10 months ago) by lgao
File size: 9047 byte(s)
Add direct solver pastix 2995. Currently works for single MPI rank only. For jobs with more than 1 MPI rank, the structure "coupleBlock" in "SystemMatrix" is not ready yet. 


1 lgao 3049
2     /*******************************************************
3     *
4     * Copyright (c) 2003-2010 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13    
14    
15     /**************************************************************/
16    
17     /* Paso: interface to the PaStiX direct library */
18    
19     /**************************************************************/
20    
21     /* Copyrights by ACcESS Australia 2010 */
22     /* Author: l.gao@uq.edu.au */
23    
24     /**************************************************************/
25    
26     #include "Paso.h"
27     #include "PASTIX.h"
28     #ifdef _OPENMP
29     #include <omp.h>
30     #endif
31    
32     /**************************************************************/
33    
34     /* free any extra stuff possibly used by the PASTIX library */
35     void Paso_PASTIX_free(Paso_SystemMatrix* A) {
36     #ifdef PASTIX
37     index_t i;
38     if (A->solver!=NULL) {
39     pastix_data_t *pastix_data = (pastix_data_t *)(A->solver);
40     pastix_int_t n = A->mainBlock->numCols;
41     pastix_int_t *perm, *invp;
42     pastix_int_t *loc2glob;
43     loc2glob = MEMALLOC(n, pastix_int_t);
44     if (A->mpi_info->size > 1)
45     for (i=0; i<n; i++) loc2glob[i]=i;
46     pastix_float_t *b;
47     pastix_int_t rhs = 1; /* number of right hand sides */
48     pastix_int_t iparm[IPARM_SIZE];
49     double dparm[DPARM_SIZE];
50     double *in;
51     perm = MEMALLOC(n, pastix_int_t);
52     invp = MEMALLOC(n, pastix_int_t);
53     in = MEMALLOC(n, double);
54    
55     /* set iparm */
56     for (i=0;i<IPARM_SIZE;++i) iparm[i]=0;
57     for (i=0;i<DPARM_SIZE;++i) dparm[i]=0.;
58     for (i=0;i<n;++i) in[i]=0.;
59    
60     iparm[IPARM_START_TASK] = API_TASK_CLEAN;
61     iparm[IPARM_END_TASK] = API_TASK_CLEAN;
62    
63     dpastix(&pastix_data, A->mpi_info->comm, n,
64     A->mainBlock->pattern->ptr, A->mainBlock->pattern->index,
65     A->mainBlock->val, loc2glob, perm, invp, in, rhs,
66     iparm, dparm);
67     MEMFREE(perm);
68     MEMFREE(invp);
69     MEMFREE(in);
70     MEMFREE(loc2glob);
71     A->solver=NULL;
72     if (iparm[IPARM_ERROR_NUMBER] != NO_ERR) Paso_setError(TYPE_ERROR,"memory release in PASTIX library failed.");
73     }
74     #endif
75     }
76    
77     /* call the solver: */
78     void Paso_PASTIX(Paso_SystemMatrix* A,
79     double* out,
80     double* in,
81     Paso_Options* options,
82     Paso_Performance* pp) {
83     #ifdef PASTIX
84     double time0, time1;
85     index_t i, m;
86     dim_t node_num, offset;
87     pastix_data_t *pastix_data = NULL;
88     pastix_int_t n = A->mainBlock->numCols;
89     pastix_int_t *perm, *invp;
90     pastix_int_t rhs = 1; /* number of right hand sides */
91     pastix_int_t iparm[IPARM_SIZE];
92     pastix_int_t *loc2glob=NULL; /* Local to global column correspondance */
93     pastix_int_t *loc2glob2=NULL; /* Local to global column correspondance */
94     pastix_int_t *ptr=NULL; /* Starting index of each column */
95     pastix_int_t *index=NULL; /* Row of each element of the matrix */
96     pastix_float_t *val=NULL; /* Value of each element of the matrix */
97     pastix_int_t *ptr2=NULL;
98     pastix_int_t *index2=NULL;
99     pastix_float_t *val2=NULL;
100     double *out2, *out_saved; /* right-hand-side */
101     double dparm[DPARM_SIZE];
102    
103     if (! (A->type & (MATRIX_FORMAT_CSC + MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
104     Paso_setError(TYPE_ERROR,"Paso_PASTIX: PASTIX requires CSC format with index offset 1 and block size 1.");
105     return;
106     }
107     options->converged=FALSE;
108     Performance_startMonitor(pp,PERFORMANCE_ALL);
109    
110     loc2glob = MEMALLOC(n,pastix_int_t);
111     offset = A->col_distribution->first_component[A->mpi_info->rank]+1;
112    
113     #pragma omp parallel for schedule(static) private(i)
114     for (i=0; i<n; i++) loc2glob[i]=i+offset;
115     /*
116     if (A->mpi_info->rank == 1){
117     fprintf(stderr, "loc2glb: ");
118     for (i=0; i<n; i++)
119     fprintf(stderr, "%d ", loc2glob[i]);
120     fprintf(stderr, "\n");
121     }
122     */
123    
124     for (i=0;i<IPARM_SIZE;++i) iparm[i]=0;
125     for (i=0;i<DPARM_SIZE;++i) dparm[i]=0.;
126    
127     perm = MEMALLOC(n, pastix_int_t);
128     invp = MEMALLOC(n, pastix_int_t);
129    
130     /* Set arrays "ptr", "index" and "val" of matrix */
131     Paso_SystemMatrix_mergeMainAndCouple(A, &ptr, &index, &val);
132    
133     /* Check Matrix format
134     * Matrix needs :
135     * - to be in fortran numbering (with offset 1)
136     * - to have only the lower triangular part in symmetric case
137     * - to have a graph with a symmetric structure in unsymmetric case
138     */
139     pastix_checkMatrix(A->mpi_info->comm, API_VERBOSE_YES,
140     ((A->type & MATRIX_FORMAT_SYM) ? API_SYM_YES:API_SYM_NO),
141     API_YES, n, &ptr, &index, &val, &loc2glob, 1);
142    
143    
144     /* initialize parameters to default values */
145     iparm[IPARM_MODIFY_PARAMETER] = API_NO;
146     dpastix(&pastix_data, A->mpi_info->comm, n,
147     ptr, index, val, loc2glob, perm, invp, in, rhs,
148     iparm, dparm);
149    
150     /* set some iparm values (customized parameters)*/
151     iparm[IPARM_VERBOSE] = API_VERBOSE_YES;
152     iparm[IPARM_ITERMAX] = 2; /* maximum number of refinements */
153     iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
154     // iparm[IPARM_RHSD_CHECK] = API_NO;
155     iparm[IPARM_START_TASK] = API_TASK_ORDERING;
156     iparm[IPARM_END_TASK] = API_TASK_BLEND;
157    
158     #ifdef _OPENMP
159     /* For thread binding on savanna, do we need to adjust the thread
160     binding to suit the structure of savanna ice node???
161     TO BE DONE */
162     /* iparm[IPARM_BINDTHRD] = API_BIND_TAB
163     pastix_setBind() */
164     iparm[IPARM_THREAD_NBR] = omp_get_max_threads();
165     #else
166     iparm[IPARM_THREAD_NBR] = 1;
167     #endif
168     iparm[IPARM_RHS_MAKING] = API_RHS_B;
169    
170     if (A->type & MATRIX_FORMAT_SYM){
171     iparm[IPARM_SYM] = API_SYM_YES;
172     iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
173     } else {
174     iparm[IPARM_SYM] = API_SYM_NO;
175     iparm[IPARM_FACTORIZATION] = API_FACT_LU;
176     }
177    
178     /* copy the right-hand-side into the solution */
179     out_saved = MEMALLOC(n, double);
180     memcpy(out_saved, in, n*sizeof(double));
181    
182     /* start solving */
183     time0=Paso_timer();
184     dpastix (&pastix_data, A->mpi_info->comm, n,
185     ptr, index, val, loc2glob, perm, invp, out_saved, 1,
186     iparm, dparm);
187    
188     node_num = pastix_getLocalNodeNbr(&pastix_data);
189     loc2glob2 = MEMALLOC(node_num, pastix_int_t);
190    
191     pastix_getLocalNodeLst(&pastix_data, loc2glob2);
192    
193     if (EXIT_SUCCESS != cscd_redispatch(n, ptr, index, val, out_saved,
194     loc2glob, node_num, &ptr2,
195     &index2, &val2, &out2, loc2glob2,
196     A->mpi_info->comm)){
197     Paso_setError(TYPE_ERROR,"Paso_PASTIX: CSCD redistribute fail.");
198     return;
199     }
200    
201     /* free part of memory allocation */
202     MEMFREE(perm);
203     MEMFREE(invp);
204     MEMFREE(ptr);
205     MEMFREE(index);
206     MEMFREE(val);
207     MEMFREE(out_saved);
208    
209     iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
210     iparm[IPARM_END_TASK] = API_TASK_CLEAN;
211    
212     dpastix(&pastix_data, A->mpi_info->comm, node_num,
213     ptr2, index2, val2, loc2glob2, perm, invp,
214     out2, 1, iparm, dparm);
215    
216     /* copy the right-hand-side back into the solution */
217     redispatch_rhs(node_num, out2, loc2glob2, n, out, loc2glob,
218     A->mpi_info->size, A->mpi_info->rank,
219     A->mpi_info->comm);
220     /*
221     {
222     MPI_Reduce(&n, &node_num, 1, MPI_INT, MPI_SUM, 0, A->mpi_info->comm);
223     if (A->mpi_info->rank == 0){
224     MEMFREE(out2);
225     out2=MEMALLOC(node_num, double);
226     }
227     MPI_Gather(out, n, MPI_DOUBLE, out2, n, MPI_DOUBLE, 0, A->mpi_info->comm);
228     if (A->mpi_info->rank == 0)
229     for (i=0; i<node_num; i++) printf("out[%d] = %g\n", i, out2[i]);
230     }*/
231    
232     options->time=Paso_timer()-time0;
233     options->set_up_time=0;
234    
235     if (iparm[IPARM_ERROR_NUMBER] != NO_ERROR) {
236     if (options->verbose) printf("PASTIX: solve failed.\n");
237     Paso_setError(VALUE_ERROR,"pastix library failed.");
238     Paso_PASTIX_free(A);
239     } else {
240     if (options->verbose) printf("PASTIX: solve completed.\n");
241     A->solver=NULL;
242     options->converged=TRUE;
243     options->residual_norm=0;
244     options->num_iter=1;
245     options->num_level=0;
246     options->num_inner_iter=0;
247     MEMFREE(loc2glob);
248     MEMFREE(loc2glob2);
249     MEMFREE(ptr2);
250     MEMFREE(index2);
251     MEMFREE(val2);
252     MEMFREE(out2);
253     }
254     Performance_stopMonitor(pp,PERFORMANCE_ALL);
255     #else
256     Paso_setError(SYSTEM_ERROR,"Paso_PASTIX: PASTIX is not avialble.");
257     #endif
258     }
259    
260     /*
261     * $Log$
262     *
263     */

  ViewVC Help
Powered by ViewVC 1.1.26