SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pzgemrdrv.c
Go to the documentation of this file.
1#include "redist.h"
2/* $Id: pzgemrdrv.c,v 1.1.1.1 2000/02/15 18:04:11 susan Exp $
3 *
4 * pzgemrdrv.c :
5 *
6 *
7 * PURPOSE:
8 *
9 * this driver is testing the PZGEMR2D routine. It calls it to obtain a new
10 * scattered block data decomposition of a distributed COMPLEX*16 (block
11 * scattered) matrix. Then it calls PZGEMR2D for the inverse redistribution
12 * and checks the results with the initial data.
13 *
14 * Data are going from a Block Scattered nbrow0 x nbcol0 decomposition on the
15 * processor grid p0 x q0, to data distributed in a BS nbrow1 x nbcol1 on the
16 * processor grid p1 x q1, then back to the BS nbrow0 x nbcol0 decomposition
17 * on the processor grid p0 x q0.
18 *
19 * See pzgemr.c file for detailed info on the PZGEMR2D function.
20 *
21 *
22 * The testing parameters are read from the file GEMR2D.dat, see the file in the
23 * distribution to have an example.
24 *
25 * created by Bernard Tourancheau in April 1994.
26 *
27 * modifications : see sccs history
28 *
29 * ===================================
30 *
31 *
32 * NOTE :
33 *
34 * - the matrix elements are COMPLEX*16
35 *
36 * - memory requirements : this procedure requires approximately 3 times the
37 * memory space of the initial data block in grid 0 (initial block, copy for
38 * test and second redistribution result) and 1 time the memory space of the
39 * result data block in grid 1. with the element size = sizeof(dcomplex)
40 * bytes,
41 *
42 *
43 * - use the procedures of the files:
44 *
45 * pzgemr.o pzgemr2.o pzgemraux.o
46 *
47 *
48 * ======================================
49 *
50 * WARNING ASSUMPTIONS :
51 *
52 *
53 * ========================================
54 *
55 *
56 * Planned changes:
57 *
58 *
59 *
60 * ========================================= */
61#define static2 static
62#if defined(Add_) || defined(f77IsF2C)
63#define fortran_mr2d pzgemr2do_
64#define fortran_mr2dnew pzgemr2d_
65#elif defined(UpCase)
66#define fortran_mr2dnew PZGEMR2D
67#define fortran_mr2d PZGEMR2DO
68#define zcopy_ ZCOPY
69#define zlacpy_ ZLACPY
70#else
71#define fortran_mr2d pzgemr2do
72#define fortran_mr2dnew pzgemr2d
73#define zcopy_ zcopy
74#define zlacpy_ zlacpy
75#endif
76#define Clacpy Czgelacpy
77void Clacpy();
78typedef struct {
79 double r, i;
80} dcomplex;
81typedef struct {
82 Int desctype;
83 Int ctxt;
84 Int m;
85 Int n;
86 Int nbrow;
87 Int nbcol;
88 Int sprow;
89 Int spcol;
90 Int lda;
91} MDESC;
92#define BLOCK_CYCLIC_2D 1
93typedef struct {
94 Int lstart;
95 Int len;
96} IDESC;
97#define SHIFT(row,sprow,nbrow) ((row)-(sprow)+ ((row) >= (sprow) ? 0 : (nbrow)))
98#define max(A,B) ((A)>(B)?(A):(B))
99#define min(A,B) ((A)>(B)?(B):(A))
100#define DIVUP(a,b) ( ((a)-1) /(b)+1)
101#define ROUNDUP(a,b) (DIVUP(a,b)*(b))
102#ifdef MALLOCDEBUG
103#define malloc mymalloc
104#define free myfree
105#define realloc myrealloc
106#endif
107/* Cblacs */
108extern void Cblacs_pcoord();
110extern void Csetpvmtids();
111extern void Cblacs_get();
112extern void Cblacs_pinfo();
113extern void Cblacs_gridinfo();
114extern void Cblacs_gridinit();
115extern void Cblacs_exit();
116extern void Cblacs_gridexit();
117extern void Cblacs_setup();
118extern void Cigebs2d();
119extern void Cigebr2d();
120extern void Cigesd2d();
121extern void Cigerv2d();
122extern void Cigsum2d();
123extern void Cigamn2d();
124extern void Cigamx2d();
125extern void Czgesd2d();
126extern void Czgerv2d();
127/* lapack */
128void zlacpy_();
129/* aux fonctions */
131extern void *mr2d_malloc();
132extern Int ppcm();
133extern Int localsize();
136extern void paramcheck();
137/* tools and others function */
138#define scanD0 zgescanD0
139#define dispmat zgedispmat
140#define setmemory zgesetmemory
141#define freememory zgefreememory
142#define scan_intervals zgescan_intervals
143extern void scanD0();
144extern void dispmat();
145extern void setmemory();
146extern void freememory();
147extern Int scan_intervals();
148extern void Cpzgemr2do();
149extern void Cpzgemr2d();
150/* some defines for Cpzgemr2do */
151#define SENDBUFF 0
152#define RECVBUFF 1
153#define SIZEBUFF 2
154#if 0
155#define DEBUG
156#endif
157#ifndef DEBUG
158#define NDEBUG
159#endif
160#include <stdio.h>
161#include <stdlib.h>
162#include <string.h>
163#include <ctype.h>
164#include <assert.h>
165/* initblock: intialize the local part of a matrix with random data (well,
166 * not very random) */
167static2 void
169{
170 dcomplex *pdata;
171 Int i;
172 pdata = block;
173 for (i = 0; i < m * n; i++, pdata++) {
174 (*pdata).r = i;
175 };
176}
177/* getparam:read from a file a list of integer parameters, the end of the
178 * parameters to read is given by a NULL at the end of the args list */
179#ifdef __STDC__
180#include <stdarg.h>
181static void
182getparam(FILE * f,...)
183{
184#else
185#include <varargs.h>
186static void
187getparam(va_alist)
188va_dcl
189{
190 FILE *f;
191#endif
192 va_list ap;
193 Int i;
194 static Int nbline;
195 char *ptr, *next;
196 Int *var;
197 static char buffer[200];
198#ifdef __STDC__
199 va_start(ap, f);
200#else
201 va_start(ap);
202 f = va_arg(ap, FILE *);
203#endif
204 do {
205 next = fgets(buffer, 200, f);
206 if (next == NULL) {
207 fprintf(stderr, "bad configuration driver file:after line %d\n", nbline);
208 exit(1);
209 }
210 nbline += 1;
211 } while (buffer[0] == '#');
212 ptr = buffer;
213 var = va_arg(ap, Int *);
214 while (var != NULL) {
215 *var = strtol(ptr, &next, 10);
216 if (ptr == next) {
217 fprintf(stderr, "bad configuration driver file:error line %d\n", nbline);
218 exit(1);
219 }
220 ptr = next;
221 var = va_arg(ap, Int *);
222 }
223 va_end(ap);
224}
225void
226initforpvm(Int argc, char *argv[])
227{
228 Int pnum, nproc;
229 Cblacs_pinfo(&pnum, &nproc);
230 if (nproc < 1) { /* we are with PVM */
231 if (pnum == 0) {
232 if (argc < 2) {
233 fprintf(stderr, "usage with PVM:xzgemr nbproc\n\
234\t where nbproc is the number of nodes to initialize\n");
235 exit(1);
236 }
237 nproc = atoi(argv[1]);
238 }
239 Cblacs_setup(&pnum, &nproc);
240 }
241}
242int
243main(int argc, char *argv[])
244{
245 /* We initialize the data-block on the current processor, then redistribute
246 * it, and perform the inverse redistribution to compare the local memory
247 * with the initial one. */
248 /* Data file */
249 FILE *fp;
250 Int nbre, nbremax;
251 /* Data distribution 0 parameters */
252 Int p0, /* # of rows in the processor grid */
253 q0; /* # of columns in the processor grid */
254 /* Data distribution 1 parameters */
255 Int p1, q1;
256 /* # of parameter to be read on the keyboard */
257#define nbparameter 24
258 /* General variables */
259 Int blocksize0;
260 Int mypnum, nprocs;
261 Int parameters[nbparameter], nberrors;
262 Int i;
263 Int ia, ja, ib, jb, m, n;
264 Int gcontext, context0, context1;
265 Int myprow1, myprow0, mypcol0, mypcol1;
266 Int dummy;
267 MDESC ma, mb;
268 dcomplex *ptrmyblock, *ptrsavemyblock, *ptrmyblockcopy, *ptrmyblockvide;
269#ifdef UsingMpiBlacs
270 MPI_Init(&argc, &argv);
271#endif
272 setvbuf(stdout, NULL, _IOLBF, 0);
273 setvbuf(stderr, NULL, _IOLBF, 0);
274#ifdef T3D
275 free(malloc(14000000));
276#endif
277 initforpvm(argc, argv);
278 /* Read physical parameters */
279 Cblacs_pinfo(&mypnum, &nprocs);
280 /* initialize BLACS for the parameter communication */
281 Cblacs_get((Int)0, (Int)0, &gcontext);
282 Cblacs_gridinit(&gcontext, "R", nprocs, (Int)1);
283 Cblacs_gridinfo(gcontext, &dummy, &dummy, &mypnum, &dummy);
284 if (mypnum == 0) {
285 if ((fp = fopen("GEMR2D.dat", "r")) == NULL) {
286 fprintf(stderr, "Can't open GEMR2D.dat\n");
287 exit(1);
288 };
289 printf("\n// ZGEMR2D TESTER for COMPLEX*16 //\n");
290 getparam(fp, &nbre, NULL);
291 printf("////////// %d tests \n\n", nbre);
292 parameters[0] = nbre;
293 Cigebs2d(gcontext, "All", "H", (Int)1, (Int)1, parameters, (Int)1);
294 } else {
295 Cigebr2d(gcontext, "All", "H", (Int)1, (Int)1, parameters, (Int)1, (Int)0, (Int)0);
296 nbre = parameters[0];
297 };
298 if (mypnum == 0) {
299 printf("\n m n m0 n0 sr0 sc0 i0 j0 p0 q0 nbr0 nbc0 \
300m1 n1 sr1 sc1 i1 j1 p1 q1 nbr1 nbc1\n\n");
301 };
302 /****** TEST LOOP *****/
303 /* Here we are in grip 1xnprocs */
304 nbremax = nbre;
305#ifdef DEBUG
306 fprintf(stderr, "bonjour,je suis le noeud %d\n", mypnum);
307#endif
308 while (nbre-- != 0) { /* Loop on the serie of tests */
309 /* All the processors read the parameters so we have to be in a 1xnprocs
310 * grid at each iteration */
311 /* Read processors grid and matrices parameters */
312 if (mypnum == 0) {
313 Int u, d;
314 getparam(fp,
315 &m, &n,
316 &ma.m, &ma.n, &ma.sprow, &ma.spcol,
317 &ia, &ja, &p0, &q0, &ma.nbrow, &ma.nbcol,
318 &mb.m, &mb.n, &mb.sprow, &mb.spcol,
319 &ib, &jb, &p1, &q1, &mb.nbrow, &mb.nbcol,
320 NULL);
321 printf("\t\t************* TEST # %d **********\n",
322 nbremax - nbre);
323 printf(" %3d %3d %3d %3d %3d %3d %3d %3d \
324%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d",
325 m, n,
326 ma.m, ma.n, ma.sprow, ma.spcol,
327 ia, ja, p0, q0, ma.nbrow, ma.nbcol,
328 mb.m, mb.n, mb.sprow, mb.spcol,
329 ib, jb, p1, q1, mb.nbrow, mb.nbcol);
330 printf("\n");
331 if (p0 * q0 > nprocs || p1 * q1 > nprocs) {
332 fprintf(stderr, "not enough nodes:%d processors required\n",
333 max(p0 * q0, p1 * q1));
334 exit(1);
335 }
336 parameters[0] = p0;
337 parameters[1] = q0;
338 parameters[2] = ma.nbrow;
339 parameters[3] = ma.nbcol;
340 parameters[4] = p1;
341 parameters[5] = q1;
342 parameters[6] = mb.nbrow;
343 parameters[7] = mb.nbcol;
344 parameters[8] = ma.m;
345 parameters[9] = ma.n;
346 parameters[10] = ma.sprow;
347 parameters[11] = ma.spcol;
348 parameters[12] = mb.sprow;
349 parameters[13] = mb.spcol;
350 parameters[14] = ia;
351 parameters[15] = ja;
352 parameters[16] = ib;
353 parameters[17] = jb;
354 parameters[18] = m;
355 parameters[19] = n;
356 parameters[20] = mb.m;
357 parameters[21] = mb.n;
358 Cigebs2d(gcontext, "All", "H", (Int)1, nbparameter, parameters, (Int)1);
359 } else {
360 Cigebr2d(gcontext, "All", "H", (Int)1, nbparameter, parameters, (Int)1, (Int)0, (Int)0);
361 p0 = parameters[0];
362 q0 = parameters[1];
363 ma.nbrow = parameters[2];
364 ma.nbcol = parameters[3];
365 p1 = parameters[4];
366 q1 = parameters[5];
367 mb.nbrow = parameters[6];
368 mb.nbcol = parameters[7];
369 ma.m = parameters[8];
370 ma.n = parameters[9];
371 ma.sprow = parameters[10];
372 ma.spcol = parameters[11];
373 mb.sprow = parameters[12];
374 mb.spcol = parameters[13];
375 ia = parameters[14];
376 ja = parameters[15];
377 ib = parameters[16];
378 jb = parameters[17];
379 m = parameters[18];
380 n = parameters[19];
381 mb.m = parameters[20];
382 mb.n = parameters[21];
385 };
386 Cblacs_get((Int)0, (Int)0, &context0);
387 Cblacs_gridinit(&context0, "R", p0, q0);
388 Cblacs_get((Int)0, (Int)0, &context1);
389 Cblacs_gridinit(&context1, "R", p1, q1);
390 Cblacs_gridinfo(context0, &dummy, &dummy, &myprow0, &mypcol0);
391 if (myprow0 >= p0 || mypcol0 >= q0)
392 myprow0 = mypcol0 = -1;
393 Cblacs_gridinfo(context1, &dummy, &dummy, &myprow1, &mypcol1);
394 if (myprow1 >= p1 || mypcol1 >= q1)
395 myprow1 = mypcol1 = -1;
396 assert((myprow0 < p0 && mypcol0 < q0) || (myprow0 == -1 && mypcol0 == -1));
397 assert((myprow1 < p1 && mypcol1 < q1) || (myprow1 == -1 && mypcol1 == -1));
398 ma.ctxt = context0;
399 mb.ctxt = context1;
400 /* From here, we are not assuming that only the processors working in the
401 * redistribution are calling xxMR2D, but the ones not concerned will do
402 * nothing. */
403 /* We compute the exact size of the local memory block for the memory
404 * allocations */
405 if (myprow0 >= 0 && mypcol0 >= 0) {
406 blocksize0 = memoryblocksize(&ma);
407 ma.lda = localsize(SHIFT(myprow0, ma.sprow, p0), p0, ma.nbrow, ma.m);
408 setmemory(&ptrmyblock, blocksize0);
409 initblock(ptrmyblock, 1, blocksize0);
410 setmemory(&ptrmyblockcopy, blocksize0);
411 memcpy((char *) ptrmyblockcopy, (char *) ptrmyblock,
412 blocksize0 * sizeof(dcomplex));
413 setmemory(&ptrmyblockvide, blocksize0);
414 for (i = 0; i < blocksize0; i++)
415 ptrmyblockvide[i].r = -1;
416 }; /* if (mypnum < p0 * q0) */
417 if (myprow1 >= 0 && mypcol1 >= 0) {
418 setmemory(&ptrsavemyblock, memoryblocksize(&mb));
419 mb.lda = localsize(SHIFT(myprow1, mb.sprow, p1), p1, mb.nbrow, mb.m);
420 }; /* if (mypnum < p1 * q1) */
421 /* Redistribute the matrix from grid 0 to grid 1 (memory location
422 * ptrmyblock to ptrsavemyblock) */
423 Cpzgemr2d(m, n,
424 ptrmyblock, ia, ja, &ma,
425 ptrsavemyblock, ib, jb, &mb, gcontext);
426 /* Perform the inverse redistribution of the matrix from grid 1 to grid 0
427 * (memory location ptrsavemyblock to ptrmyblockvide) */
428 Cpzgemr2d(m, n,
429 ptrsavemyblock, ib, jb, &mb,
430 ptrmyblockvide, ia, ja, &ma, gcontext);
431 /* Check the differences */
432 nberrors = 0;
433 if (myprow0 >= 0 && mypcol0 >= 0) {
434 /* only for the processors that do have data at the begining */
435 for (i = 0; i < blocksize0; i++) {
436 Int li, lj, gi, gj;
437 Int in;
438 in = 1;
439 li = i % ma.lda;
440 lj = i / ma.lda;
441 gi = (li / ma.nbrow) * p0 * ma.nbrow +
442 SHIFT(myprow0, ma.sprow, p0) * ma.nbrow + li % ma.nbrow;
443 gj = (lj / ma.nbcol) * q0 * ma.nbcol +
444 SHIFT(mypcol0, ma.spcol, q0) * ma.nbcol + lj % ma.nbcol;
445 assert(gi < ma.m && gj < ma.n);
446 gi -= (ia - 1);
447 gj -= (ja - 1);
448 if (gi < 0 || gj < 0 || gi >= m || gj >= n)
449 in = 0;
450 if (!in) {
451 ptrmyblockcopy[i].r = -1;
452 }
453 if (ptrmyblockvide[i].r != ptrmyblockcopy[i].r) {
454 nberrors++;
455 printf("Proc %d : Error element number %d, value = %f , initvalue =%f \n"
456 ,mypnum, i,
457 ptrmyblockvide[i].r, ptrmyblockcopy[i].r);
458 };
459 };
460 if (nberrors > 0) {
461 printf("Processor %d, has tested %d COMPLEX*16 elements,\
462Number of redistribution errors = %d \n",
463 mypnum, blocksize0, nberrors);
464 }
465 }
466 /* Look at the errors on all the processors at this point. */
467 Cigsum2d(gcontext, "All", "H", (Int)1, (Int)1, &nberrors, (Int)1, (Int)0, (Int)0);
468 if (mypnum == 0)
469 if (nberrors)
470 printf(" => Total number of redistribution errors = %d \n",
471 nberrors);
472 else
473 printf("TEST PASSED OK\n");
474 /* release memory for the next iteration */
475 if (myprow0 >= 0 && mypcol0 >= 0) {
476 freememory((char *) ptrmyblock);
477 freememory((char *) ptrmyblockvide);
478 freememory((char *) ptrmyblockcopy);
479 }; /* if (mypnum < p0 * q0) */
480 /* release memory for the next iteration */
481 if (myprow1 >= 0 && mypcol1 >= 0) {
482 freememory((char *) ptrsavemyblock);
483 };
484 if (myprow0 >= 0)
485 Cblacs_gridexit(context0);
486 if (myprow1 >= 0)
487 Cblacs_gridexit(context1);
488 }; /* while nbre != 0 */
489 if (mypnum == 0) {
490 fclose(fp);
491 };
492 Cblacs_exit((Int)0);
493 return 0;
494}/* main */
#define Int
Definition Bconfig.h:22
Int memoryblocksize()
Int changeorigin()
#define freememory
Definition pzgemrdrv.c:141
void Czgesd2d()
#define scan_intervals
Definition pzgemrdrv.c:142
#define SHIFT(row, sprow, nbrow)
Definition pzgemrdrv.c:97
void Cblacs_gridexit()
#define max(A, B)
Definition pzgemrdrv.c:98
void Cigsum2d()
#define static2
Definition pzgemrdrv.c:61
void Cpzgemr2d()
#define scanD0
Definition pzgemrdrv.c:138
Int Cblacs_pnum()
void Cpzgemr2do()
#define Clacpy
Definition pzgemrdrv.c:76
void Cigamx2d()
void Cblacs_pinfo()
void Cblacs_pcoord()
void Cigamn2d()
Int localsize()
#define setmemory
Definition pzgemrdrv.c:140
void Cblacs_get()
#define zlacpy_
Definition pzgemrdrv.c:74
void Czgerv2d()
void Cigerv2d()
void Cigebs2d()
#define dispmat
Definition pzgemrdrv.c:139
static2 void initblock(dcomplex *block, Int m, Int n)
Definition pzgemrdrv.c:168
void paramcheck()
void Cblacs_setup()
#define BLOCK_CYCLIC_2D
Definition pzgemrdrv.c:92
void Cblacs_gridinit()
void Cblacs_gridinfo()
void Csetpvmtids()
Int localindice()
void Cigesd2d()
void initforpvm(Int argc, char *argv[])
Definition pzgemrdrv.c:226
Int ppcm()
void Cigebr2d()
void Cblacs_exit()
void * mr2d_malloc()
#define nbparameter
main()
Definition size.c:2
Int m
Definition pcgemr.c:166
Int spcol
Definition pcgemr.c:171
Int nbcol
Definition pcgemr.c:169
Int sprow
Definition pcgemr.c:170
Int nbrow
Definition pcgemr.c:168
Int ctxt
Definition pcgemr.c:165
Int desctype
Definition pcgemr.c:164
Int n
Definition pcgemr.c:167
Int lda
Definition pcgemr.c:172
double r
Definition pzgemr.c:161