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

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

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 MIME type: text/plain
File size: 8591 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 iparm[IPARM_SIZE];
48 double dparm[DPARM_SIZE];
49 double *in;
50
51 perm = MEMALLOC(n, pastix_int_t);
52 invp = MEMALLOC(n, pastix_int_t);
53 in = MEMALLOC(n, double);
54
55 for (i=0;i<IPARM_SIZE;++i) iparm[i]=0;
56 for (i=0;i<DPARM_SIZE;++i) dparm[i]=0.;
57 for (i=0;i<n;++i) in[i]=0.;
58
59 iparm[IPARM_START_TASK] = API_TASK_CLEAN;
60 iparm[IPARM_END_TASK] = API_TASK_CLEAN;
61
62 dpastix(&pastix_data, A->mpi_info->comm, n,
63 A->mainBlock->pattern->ptr, A->mainBlock->pattern->index,
64 A->mainBlock->val, loc2glob, perm, invp, in, 1,
65 iparm, dparm);
66 MEMFREE(perm);
67 MEMFREE(invp);
68 MEMFREE(in);
69 MEMFREE(loc2glob);
70 A->solver=NULL;
71 if (iparm[IPARM_ERROR_NUMBER] != NO_ERR) Paso_setError(TYPE_ERROR,"memory release in PASTIX library failed.");
72 }
73 #endif
74 }
75
76 /* call the solver: */
77 void Paso_PASTIX(Paso_SystemMatrix* A,
78 double* out,
79 double* in,
80 Paso_Options* options,
81 Paso_Performance* pp) {
82 #ifdef PASTIX
83 double time0;
84 index_t i, m;
85 dim_t node_num, offset;
86 pastix_data_t *pastix_data = NULL;
87 pastix_int_t n=A->mainBlock->numCols;
88 pastix_int_t *perm=NULL;
89 pastix_int_t *invp=NULL;
90 pastix_int_t iparm[IPARM_SIZE];
91 pastix_int_t *loc2glob=NULL; /* Local to global column correspondance */
92 pastix_int_t *loc2glob2=NULL; /* Local to global column correspondance */
93 pastix_int_t *ptr=NULL; /* Starting index of each column */
94 pastix_int_t *index=NULL; /* Row of each element of the matrix */
95 pastix_float_t *val=NULL; /* Value of each element of the matrix */
96 pastix_int_t *ptr2=NULL;
97 pastix_int_t *index2=NULL;
98 pastix_float_t *val2=NULL;
99 double *out_saved=NULL; /* right-hand-side */
100 double *out2=NULL;
101 double dparm[DPARM_SIZE];
102
103 printf("CHP 0\n");
104
105 if (! (A->type & (MATRIX_FORMAT_CSC + MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) ) {
106 Paso_setError(TYPE_ERROR,"Paso_PASTIX: PASTIX requires CSC format with index offset 1 and block size 1.");
107 return;
108 }
109 options->converged=FALSE;
110 Performance_startMonitor(pp,PERFORMANCE_ALL);
111
112 loc2glob = MEMALLOC(n,pastix_int_t);
113 perm = MEMALLOC(n, pastix_int_t);
114 invp = MEMALLOC(n, pastix_int_t);
115
116 offset = A->col_distribution->first_component[A->mpi_info->rank]+1;
117 #pragma omp parallel for schedule(static) private(i)
118 for (i=0; i<n; i++) loc2glob[i]=i+offset;
119 #pragma omp parallel for schedule(static) private(i)
120 for (i=0;i<IPARM_SIZE;++i) iparm[i]=0;
121 #pragma omp parallel for schedule(static) private(i)
122 for (i=0;i<DPARM_SIZE;++i) dparm[i]=0.;
123
124 /* Set arrays "ptr", "index" and "val" of matrix */
125 Paso_SystemMatrix_mergeMainAndCouple(A, &ptr, &index, &val);
126
127 /* Check Matrix format
128 * Matrix needs :
129 * - to be in fortran numbering (with offset 1)
130 * - to have only the lower triangular part in symmetric case
131 * - to have a graph with a symmetric structure in unsymmetric case
132 */
133 pastix_checkMatrix(A->mpi_info->comm, API_VERBOSE_YES,
134 ((A->type & MATRIX_FORMAT_SYM) ? API_SYM_YES:API_SYM_NO),
135 API_YES, n, &ptr, &index, &val, &loc2glob, 1);
136
137 /* copy the right-hand-side into the solution */
138 out_saved = MEMALLOC(n, double);
139 memcpy(out_saved, in, n*sizeof(double));
140
141 /* initialize parameters to default values */
142 iparm[IPARM_MODIFY_PARAMETER] = API_NO;
143 dpastix(&pastix_data, A->mpi_info->comm, n,
144 ptr, index, val, loc2glob, perm, invp, out_saved, 1,
145 iparm, dparm);
146
147 /* set some iparm values (customized parameters)*/
148 iparm[IPARM_VERBOSE] = API_VERBOSE_YES;
149 iparm[IPARM_ITERMAX] = 2; /* maximum number of refinements */
150 iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
151 /* iparm[IPARM_RHSD_CHECK] = API_NO; */
152 iparm[IPARM_START_TASK] = API_TASK_ORDERING;
153 iparm[IPARM_END_TASK] = API_TASK_BLEND;
154 #ifdef _OPENMP
155 /* For thread binding on savanna, do we need to adjust the thread
156 binding to suit the structure of savanna ice node???
157 TO BE DONE */
158 /* iparm[IPARM_BINDTHRD] = API_BIND_TAB
159 pastix_setBind() */
160 iparm[IPARM_THREAD_NBR] = omp_get_max_threads();
161 #else
162 iparm[IPARM_THREAD_NBR] = 1;
163 #endif
164 iparm[IPARM_RHS_MAKING] = API_RHS_B;
165 if (A->type & MATRIX_FORMAT_SYM){
166 iparm[IPARM_SYM] = API_SYM_YES;
167 iparm[IPARM_FACTORIZATION] = API_FACT_LDLT;
168 } else {
169 iparm[IPARM_SYM] = API_SYM_NO;
170 iparm[IPARM_FACTORIZATION] = API_FACT_LU;
171 }
172
173 /* start solving */
174 time0=Paso_timer();
175 dpastix (&pastix_data, A->mpi_info->comm, n,
176 ptr, index, val, loc2glob, perm, invp, out_saved, 1,
177 iparm, dparm);
178
179 node_num = pastix_getLocalNodeNbr(&pastix_data);
180 loc2glob2 = MEMALLOC(node_num, pastix_int_t);
181
182 pastix_getLocalNodeLst(&pastix_data, loc2glob2);
183
184 if (EXIT_SUCCESS != cscd_redispatch(n, ptr, index, val, out_saved,
185 loc2glob, node_num, &ptr2,
186 &index2, &val2, &out2, loc2glob2,
187 A->mpi_info->comm)){
188 Paso_setError(TYPE_ERROR,"Paso_PASTIX: CSCD redistribute fail.");
189 return;
190 }
191
192 /* free part of memory allocation */
193 MEMFREE(perm);
194 MEMFREE(invp);
195 MEMFREE(ptr);
196 MEMFREE(index);
197 MEMFREE(val);
198 MEMFREE(out_saved);
199
200 iparm[IPARM_START_TASK] = API_TASK_NUMFACT;
201 iparm[IPARM_END_TASK] = API_TASK_CLEAN;
202
203 dpastix(&pastix_data, A->mpi_info->comm, node_num,
204 ptr2, index2, val2, loc2glob2, perm, invp,
205 out2, 1, iparm, dparm);
206
207 /* copy the right-hand-side back into the solution */
208 redispatch_rhs(node_num, out2, loc2glob2, n, out, loc2glob,
209 A->mpi_info->size, A->mpi_info->rank,
210 A->mpi_info->comm);
211
212 options->time=Paso_timer()-time0;
213 options->set_up_time=0;
214
215 if (iparm[IPARM_ERROR_NUMBER] != NO_ERROR) {
216 if (options->verbose) printf("PASTIX: solve failed.\n");
217 Paso_setError(VALUE_ERROR,"pastix library failed.");
218 Paso_PASTIX_free(A);
219 } else {
220 if (options->verbose) printf("PASTIX: solve completed.\n");
221 // A->solver=NULL;
222 options->converged=TRUE;
223 options->residual_norm=0;
224 options->num_iter=1;
225 options->num_level=0;
226 options->num_inner_iter=0;
227 MEMFREE(loc2glob);
228 MEMFREE(loc2glob2);
229 MEMFREE(ptr2);
230 MEMFREE(index2);
231 MEMFREE(val2);
232 MEMFREE(out2);
233 }
234 Performance_stopMonitor(pp,PERFORMANCE_ALL);
235 #else
236 Paso_setError(SYSTEM_ERROR,"Paso_PASTIX: PASTIX is not avialble.");
237 #endif
238 }
239
240 /*
241 * $Log$
242 *
243 */

  ViewVC Help
Powered by ViewVC 1.1.26