1 /* testMPI.c */
2
3 #include "../BridgeMPI.h"
4
5 void JimMatMulMPI ( ) ;
6 void JimSolveMPI ( ) ;
7
8 /*--------------------------------------------------------------------*/
9
main(int argc,char * argv[])10 void main ( int argc, char *argv[] )
11 /*
12 -----------------------------------------------------------------
13 MPI environment: read in Harwell-Boeing matrices, using factor,
14 solve, and multiply routines based on spooles, invoke eigensolver
15
16 created -- 98mar31, jcp
17 modified -- 98dec18, cca
18 -----------------------------------------------------------------
19 */
20 {
21 BridgeMPI bridge ;
22 MPI_Comm comm ;
23
24 char *inFileName_A, *inFileName_B, *parmFileName, *type ;
25 char buffer[20], pbtype[4], which[4] ;
26 int error, fstevl, lfinit, lstevl, msglvl, myid, mxbksz, ncol,
27 ndiscd, neig, neigvl, nfound, nnonzeros, nproc, nrhs, nrow,
28 prbtyp, rc, retc, rfinit, seed, warnng ;
29 int c__5 = 5, output = 6 ;
30 int *lanczos_wksp ;
31 InpMtx *inpmtxA, *inpmtxB ;
32 FILE *msgFile, *parmFile ;
33 double lftend, rhtend, center, shfscl, t1, t2 ;
34 double c__1 = 1.0, c__4 = 4.0, tolact = 2.309970868130169e-11 ;
35 double eigval[1000], sigma[2] ;
36 double *evec;
37 /*
38 ---------------------------------------------------------------
39 find out the identity of this process and the number of process
40 ---------------------------------------------------------------
41 */
42 MPI_Init(&argc, &argv) ;
43 MPI_Comm_dup(MPI_COMM_WORLD, &comm) ;
44 MPI_Comm_rank(comm, &myid) ;
45 MPI_Comm_size(comm, &nproc) ;
46 fprintf(stdout, "\n myid = %d", myid) ;
47 fflush(stdout) ;
48 /*--------------------------------------------------------------------*/
49 /*
50 -----------------------------
51 decode the command line input
52 -----------------------------
53 */
54 if ( argc != 7 ) {
55 fprintf(stdout,
56 "\n\n usage : %s msglvl msgFile parmFile seed inFileA inFileB"
57 "\n msglvl -- message level"
58 "\n msgFile -- message file"
59 "\n parmFile -- input parameters file"
60 "\n seed -- random number seed, used for ordering"
61 "\n inFileA -- stiffness matrix, in Harwell-Boeing format"
62 "\n inFileB -- mass matrix, in Harwell-Boeing format"
63 "\n used for prbtyp = 1 or 2"
64 "\n", argv[0]) ;
65 return ;
66 }
67 msglvl = atoi(argv[1]) ;
68 if ( strcmp(argv[2], "stdout") == 0 ) {
69 msgFile = stdout ;
70 } else {
71 int length = strlen(argv[2]) + 1 + 4 ;
72 char *buffer = CVinit(length, '\0') ;
73 sprintf(buffer, "%s.%d", argv[2], myid) ;
74 if ( (msgFile = fopen(buffer, "w")) == NULL ) {
75 fprintf(stderr, "\n fatal error in %s"
76 "\n unable to open file %s\n",
77 argv[0], buffer) ;
78 return ;
79 }
80 CVfree(buffer) ;
81 }
82 parmFileName = argv[3] ;
83 seed = atoi(argv[4]) ;
84 inFileName_A = argv[5] ;
85 inFileName_B = argv[6] ;
86 fprintf(msgFile,
87 "\n %s "
88 "\n msglvl -- %d"
89 "\n message file -- %s"
90 "\n parameter file -- %s"
91 "\n stiffness matrix file -- %s"
92 "\n mass matrix file -- %s"
93 "\n random number seed -- %d"
94 "\n",
95 argv[0], msglvl, argv[2], parmFileName, inFileName_A,
96 inFileName_B, seed) ;
97 fflush(msgFile) ;
98 if ( strcmp(inFileName_A, "none") == 0 ) {
99 fprintf(msgFile, "\n no file to read from") ;
100 exit(0) ;
101 }
102 /*--------------------------------------------------------------------*/
103 if ( myid == 0 ) {
104 /*
105 ----------------------------------------------
106 processor zero reads in the matrix header info
107 ----------------------------------------------
108 */
109 MARKTIME(t1) ;
110 readHB_info(inFileName_A, &nrow, &ncol, &nnonzeros, &type, &nrhs) ;
111 MARKTIME(t2) ;
112 fprintf(msgFile, "\n CPU %8.3f : read in harwell-boeing header info",
113 t2 - t1) ;
114 fflush(msgFile) ;
115 }
116 MPI_Bcast((void *) &nrow, 1, MPI_INT, 0, MPI_COMM_WORLD) ;
117 /*--------------------------------------------------------------------*/
118 /*
119 ---------------------------------------------------------------
120 read in eigenvalue problem data
121 neigvl -- # of desired eigenvalues
122 which -- which eigenvalues to compute
123 'l' or 'L' lowest (smallest magnitude)
124 'h' or 'H' highest (largest magnitude)
125 'n' or 'N' nearest to central value
126 'c' or 'C' nearest to central value
127 'a' or 'A' all eigenvalues in interval
128 pbtype -- type of problem
129 'v' or 'V' generalized symmetric problem (K,M)
130 with M positive semidefinite (vibration problem)
131 'b' or 'B' generalized symmetric problem (K,K_s)
132 with K positive semidefinite
133 with K_s posibly indefinite (buckling problem)
134 'o' or 'O' ordinary symmetric eigenproblem
135 lfinit -- if true, lftend is restriction on lower bound of
136 eigenvalues. if false, no restriction on lower bound
137 lftend -- left endpoint of interval
138 rfinit -- if true, rhtend is restriction on upper bound of
139 eigenvalues. if false, no restriction on upper bound
140 rhtend -- right endpoint of interval
141 center -- center of interval
142 mxbksz -- upper bound on block size for Lanczos recurrence
143 shfscl -- shift scaling parameter, an estimate on the magnitude
144 of the smallest nonzero eigenvalues
145 ---------------------------------------------------------------
146 */
147 MARKTIME(t1) ;
148 parmFile = fopen(parmFileName, "r");
149 fscanf(parmFile, "%d %s %s %d %le %d %le %le %d %le",
150 &neigvl, which, pbtype, &lfinit, &lftend,
151 &rfinit, &rhtend, ¢er, &mxbksz, &shfscl) ;
152 fclose(parmFile);
153 MARKTIME(t2) ;
154 fprintf(msgFile, "\n CPU %8.3f : read in eigenvalue problem data",
155 t2 - t1) ;
156 fflush(msgFile) ;
157 /*
158 ----------------------------------------
159 check and set the problem type parameter
160 ----------------------------------------
161 */
162 switch ( pbtype[1] ) {
163 case 'v' : case 'V' : prbtyp = 1 ; break ;
164 case 'b' : case 'B' : prbtyp = 2 ; break ;
165 case 'o' : case 'O' : prbtyp = 3 ; break ;
166 default :
167 fprintf(stderr, "\n invalid problem type %s", pbtype) ;
168 exit(-1) ;
169 }
170 /*
171 ----------------------------
172 Initialize Lanczos workspace
173 ----------------------------
174 */
175 MARKTIME(t1) ;
176 lanczos_init_ ( &lanczos_wksp ) ;
177 MARKTIME(t2) ;
178 fprintf(msgFile, "\n CPU %8.3f : initialize Lanczos workspace",
179 t2 - t1) ;
180 fflush(msgFile) ;
181 /*
182 ----------------------------------
183 initialize communication structure
184 ----------------------------------
185 */
186 MARKTIME(t1) ;
187 lanczos_set_parm( &lanczos_wksp, "order-of-problem", &nrow, &retc );
188 lanczos_set_parm( &lanczos_wksp, "accuracy-tolerance", &tolact, &retc );
189 lanczos_set_parm( &lanczos_wksp, "max-block-size", &mxbksz, &retc );
190 lanczos_set_parm( &lanczos_wksp, "shift-scale", &shfscl, &retc );
191 lanczos_set_parm( &lanczos_wksp, "message_level", &msglvl, &retc );
192 lanczos_set_parm( &lanczos_wksp, "mpi-communicator", &comm, &retc );
193 lanczos_set_parm( &lanczos_wksp, "qfile-pathname", "lqfil", &retc );
194 lanczos_set_parm( &lanczos_wksp, "mqfil-pathname", "lmqfil", &retc );
195 lanczos_set_parm( &lanczos_wksp, "evfil-pathname", "evcfil", &retc );
196 MARKTIME(t2) ;
197 fprintf(msgFile,
198 "\n CPU %8.3f : init the lanczos communication structure",
199 t2 - t1) ;
200 fflush(msgFile) ;
201 /*--------------------------------------------------------------------*/
202 if ( myid == 0 ) {
203 /*
204 ------------------------------------
205 processor zero reads in the matrices
206 ------------------------------------
207 */
208 MARKTIME(t1) ;
209 inpmtxA = InpMtx_new() ;
210 InpMtx_readFromHBfile ( inpmtxA, inFileName_A ) ;
211 MARKTIME(t2) ;
212 fprintf(msgFile, "\n CPU %8.3f : read in first matrix", t2 - t1) ;
213 fflush(msgFile) ;
214 if ( msglvl > 2 ) {
215 fprintf(msgFile, "\n\n InpMtx A object after loading") ;
216 InpMtx_writeForHumanEye(inpmtxA, msgFile) ;
217 fflush(msgFile) ;
218 }
219 lanczos_set_parm( &lanczos_wksp, "matrix-type", &c__1, &retc );
220 if ( prbtyp != 3 ) {
221 if ( strcmp(inFileName_B, "none") == 0 ) {
222 fprintf(msgFile, "\n no file to read from") ;
223 exit(0) ;
224 }
225 MARKTIME(t1) ;
226 inpmtxB = InpMtx_new() ;
227 InpMtx_readFromHBfile ( inpmtxB, inFileName_B ) ;
228 MARKTIME(t2) ;
229 fprintf(msgFile, "\n CPU %8.3f : read in first matrix", t2 - t1) ;
230 fflush(msgFile) ;
231 if ( msglvl > 2 ) {
232 fprintf(msgFile, "\n\n InpMtx B object after loading") ;
233 InpMtx_writeForHumanEye(inpmtxB, msgFile) ;
234 fflush(msgFile) ;
235 }
236 } else {
237 inpmtxB = NULL ;
238 lanczos_set_parm( &lanczos_wksp, "matrix-type", &c__4, &retc );
239 }
240 } else {
241 /*
242 ------------------------------------------------
243 other processors initialize their local matrices
244 ------------------------------------------------
245 */
246 inpmtxA = InpMtx_new() ;
247 InpMtx_init(inpmtxA, INPMTX_BY_CHEVRONS, SPOOLES_REAL, 0, 0) ;
248 lanczos_set_parm( &lanczos_wksp, "matrix-type", &c__1, &retc );
249 if ( prbtyp == 1 || prbtyp == 2 ) {
250 inpmtxB = InpMtx_new() ;
251 InpMtx_init(inpmtxB, INPMTX_BY_CHEVRONS, SPOOLES_REAL, 0, 0) ;
252 } else {
253 inpmtxB = NULL ;
254 lanczos_set_parm( &lanczos_wksp, "matrix-type", &c__4, &retc );
255 }
256 }
257 /*
258 -----------------------------
259 set up the solver environment
260 -----------------------------
261 */
262 MARKTIME(t1) ;
263 rc = SetupMPI((void *) &bridge, &prbtyp, &nrow, &mxbksz, inpmtxA,
264 inpmtxB, &seed, &msglvl, msgFile, MPI_COMM_WORLD) ;
265 MARKTIME(t2) ;
266 fprintf(msgFile, "\n CPU %8.3f : set up solver environment", t2 - t1) ;
267 fflush(msgFile) ;
268 if ( rc != 1 ) {
269 fprintf(stderr, "\n fatal error return %d from SetupMPI()", rc) ;
270 MPI_Finalize() ;
271 exit(-1) ;
272 }
273 /*--------------------------------------------------------------------*/
274 /*
275 -----------------------------------------------
276 invoke eigensolver
277 nfound -- # of eigenvalues found and kept
278 ndisc -- # of additional eigenvalues discarded
279 -----------------------------------------------
280 */
281 MARKTIME(t1) ;
282 lanczos_run ( &neigvl, &which[1] , &pbtype[1], &lfinit, &lftend,
283 &rfinit, &rhtend, ¢er, &lanczos_wksp, &bridge, &nfound,
284 &ndiscd, &warnng, &error, FactorMPI, JimMatMulMPI,
285 JimSolveMPI ) ;
286 MARKTIME(t2) ;
287 fprintf(msgFile, "\n CPU %8.3f : time for lanczos run", t2 - t1) ;
288 fflush(msgFile) ;
289 if ( myid == 0 ) {
290 /*
291 ----------------------------------------------
292 processor 0 deals with eigenvalues and vectors
293 ----------------------------------------------
294 */
295 MARKTIME(t1) ;
296 neig = nfound + ndiscd ;
297 lstevl = nfound ;
298 lanczos_eigenvalues (&lanczos_wksp, eigval, &neig, &retc);
299 fstevl = 1 ;
300 if ( nfound == 0 ) fstevl = -1 ;
301 if ( ndiscd > 0 ) lstevl = -ndiscd ;
302 hdslp5_ ("computed eigenvalues returned by hdserl",
303 &neig, eigval, &output, 39L ) ;
304 MARKTIME(t2) ;
305 fprintf(msgFile, "\n CPU %8.3f : get and print eigenvalues",
306 t2 - t1) ;
307 fflush(msgFile) ;
308 /*
309 -------------------------
310 get eigenvectors and print
311 -------------------------
312 */
313 /*
314 MARKTIME(t1) ;
315 neig = min ( 50, nrow );
316 Lncz_ALLOCATE(evec, double, nrow, retc);
317 for (i = 1; i<= nfound; i++) {
318 lanczos_eigenvector(&lanczos_wksp, &i, &i, newToOld,
319 evec, &nrow, &retc) ;
320 hdslp5_("computed eigenvector returned by hdserc",
321 &neig, evec, &output, 39L ) ;
322 }
323 MARKTIME(t2) ;
324 fprintf(msgFile, "\n CPU %8.3f : get and print eigenvectors",
325 t2 - t1) ;
326 fflush(msgFile) ;
327 */
328 }
329 /*
330 ------------------------
331 free the working storage
332 ------------------------
333 */
334 MARKTIME(t1) ;
335 lanczos_free(&lanczos_wksp) ;
336 MARKTIME(t2) ;
337 fprintf(msgFile, "\n CPU %8.3f : free lanczos workspace", t2 - t1) ;
338 fflush(msgFile) ;
339 MARKTIME(t1) ;
340 CleanupMPI(&bridge) ;
341 MARKTIME(t2) ;
342 fprintf(msgFile, "\n CPU %8.3f : free solver workspace", t2 - t1) ;
343 fflush(msgFile) ;
344
345 MPI_Finalize() ;
346
347 fprintf(msgFile, "\n") ;
348 fclose(msgFile) ;
349
350 return ; }
351
352