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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3049 - (show 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
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