SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pstrmr.c
Go to the documentation of this file.
1#include "redist.h"
158#define static2 static
159#if defined(Add_) || defined(f77IsF2C)
160#define fortran_mr2d pstrmr2do_
161#define fortran_mr2dnew pstrmr2d_
162#elif defined(UpCase)
163#define fortran_mr2dnew PSTRMR2D
164#define fortran_mr2d PSTRMR2DO
165#define scopy_ SCOPY
166#define slacpy_ SLACPY
167#else
168#define fortran_mr2d pstrmr2do
169#define fortran_mr2dnew pstrmr2d
170#define scopy_ scopy
171#define slacpy_ slacpy
172#endif
173#define Clacpy Cstrlacpy
174void Clacpy();
175typedef struct {
176 Int desctype;
177 Int ctxt;
178 Int m;
179 Int n;
180 Int nbrow;
181 Int nbcol;
182 Int sprow;
183 Int spcol;
184 Int lda;
185} MDESC;
186#define BLOCK_CYCLIC_2D 1
187typedef struct {
188 Int gstart;
189 Int len;
190} IDESC;
191#define SHIFT(row,sprow,nbrow) ((row)-(sprow)+ ((row) >= (sprow) ? 0 : (nbrow)))
192#define max(A,B) ((A)>(B)?(A):(B))
193#define min(A,B) ((A)>(B)?(B):(A))
194#define DIVUP(a,b) ( ((a)-1) /(b)+1)
195#define ROUNDUP(a,b) (DIVUP(a,b)*(b))
196#ifdef MALLOCDEBUG
197#define malloc mymalloc
198#define free myfree
199#define realloc myrealloc
200#endif
201/* Cblacs */
202extern void Cblacs_pcoord();
204extern void Csetpvmtids();
205extern void Cblacs_get();
206extern void Cblacs_pinfo();
207extern void Cblacs_gridinfo();
208extern void Cblacs_gridinit();
209extern void Cblacs_exit();
210extern void Cblacs_gridexit();
211extern void Cblacs_setup();
212extern void Cigebs2d();
213extern void Cigebr2d();
214extern void Cigesd2d();
215extern void Cigerv2d();
216extern void Cigsum2d();
217extern void Cigamn2d();
218extern void Cigamx2d();
219extern void Csgesd2d();
220extern void Csgerv2d();
221/* lapack */
222void slacpy_();
223/* aux fonctions */
225extern void *mr2d_malloc();
226extern Int ppcm();
227extern Int localsize();
230extern void paramcheck();
231/* tools and others function */
232#define scanD0 strscanD0
233#define dispmat strdispmat
234#define setmemory strsetmemory
235#define freememory strfreememory
236#define scan_intervals strscan_intervals
237extern void scanD0();
238extern void dispmat();
239extern void setmemory();
240extern void freememory();
241extern Int scan_intervals();
242extern void Cpstrmr2do();
243extern void Cpstrmr2d();
244/* some defines for Cpstrmr2do */
245#define SENDBUFF 0
246#define RECVBUFF 1
247#define SIZEBUFF 2
248#if 0
249#define DEBUG
250#endif
251#ifndef DEBUG
252#define NDEBUG
253#endif
254#include <stdio.h>
255#include <stdlib.h>
256#include <assert.h>
257#define DESCLEN 9
258void
259fortran_mr2d(char *uplo, char *diag, Int *m, Int *n, float *A, Int *ia, Int *ja, Int desc_A[DESCLEN],
260 float *B, Int *ib, Int *jb, Int desc_B[DESCLEN])
261{
262 Cpstrmr2do(uplo, diag, *m, *n, A, *ia, *ja, (MDESC *) desc_A,
263 B, *ib, *jb, (MDESC *) desc_B);
264 return;
265}
266void
267fortran_mr2dnew(char *uplo, char *diag, Int *m, Int *n, float *A, Int *ia, Int *ja, Int desc_A[DESCLEN],
268 float *B, Int *ib, Int *jb, Int desc_B[DESCLEN], Int *gcontext)
269{
270 Cpstrmr2d(uplo, diag, *m, *n, A, *ia, *ja, (MDESC *) desc_A,
271 B, *ib, *jb, (MDESC *) desc_B, *gcontext);
272 return;
273}
279void
280Cpstrmr2do(uplo, diag, m, n,
281 ptrmyblock, ia, ja, ma,
282 ptrmynewblock, ib, jb, mb)
283 char *uplo, *diag;
284 float *ptrmyblock, *ptrmynewblock;
285/* pointers to the memory location of the matrix and the redistributed matrix */
286 MDESC *ma;
287 MDESC *mb;
288 Int ia, ja, ib, jb, m, n;
289{
290 Int dummy, nprocs;
291 Int gcontext;
292 /* first we initialize a global grid which serve as a reference to
293 * communicate from grid a to grid b */
294 Cblacs_pinfo(&dummy, &nprocs);
295 Cblacs_get((Int)0, (Int)0, &gcontext);
296 Cblacs_gridinit(&gcontext, "R", (Int)1, nprocs);
297 Cpstrmr2d(uplo, diag, m, n, ptrmyblock, ia, ja, ma,
298 ptrmynewblock, ib, jb, mb, gcontext);
299 Cblacs_gridexit(gcontext);
300}
301#define NBPARAM 20 /* p0,q0,p1,q1, puis ma,na,mba,nba,rowa,cola puis
302 * idem B puis ia,ja puis ib,jb */
303#define MAGIC_MAX 100000000
304void
305Cpstrmr2d(char *uplo, char *diag, Int m, Int n,
306 float *ptrmyblock, Int ia, Int ja, MDESC *ma,
307 float *ptrmynewblock, Int ib, Int jb, MDESC *mb, Int globcontext)
308{
309 float *ptrsendbuff, *ptrrecvbuff, *ptrNULL = 0;
310 float *recvptr;
311 MDESC newa, newb;
312 Int *proc0, *proc1, *param;
313 Int mypnum, myprow0, mypcol0, myprow1, mypcol1, nprocs;
314 Int i, j;
315 Int nprow, npcol, gcontext;
316 Int recvsize, sendsize;
317 IDESC *h_inter; /* to store the horizontal intersections */
318 IDESC *v_inter; /* to store the vertical intersections */
319 Int hinter_nb, vinter_nb; /* number of intrsections in both directions */
320 Int dummy;
321 Int p0, q0, p1, q1;
322 Int *ra, *ca;
323 /* end of variables */
324 /* To simplify further calcul we change the matrix indexation from
325 * 1..m,1..n (fortran) to 0..m-1,0..n-1 */
326 if (m == 0 || n == 0)
327 return;
328 ia -= 1;
329 ja -= 1;
330 ib -= 1;
331 jb -= 1;
332 Cblacs_gridinfo(globcontext, &nprow, &npcol, &dummy, &mypnum);
333 gcontext = globcontext;
334 nprocs = nprow * npcol;
335 /* if the global context that is given to us has not the shape of a line
336 * (nprow != 1), create a new context. TODO: to be optimal, we should
337 * avoid this because it is an uncessary synchronisation */
338 if (nprow != 1) {
339 gridreshape(&gcontext);
340 Cblacs_gridinfo(gcontext, &dummy, &dummy, &dummy, &mypnum);
341 }
342 Cblacs_gridinfo(ma->ctxt, &p0, &q0, &myprow0, &mypcol0);
343 /* compatibility T3D, must check myprow and mypcol are within bounds */
344 if (myprow0 >= p0 || mypcol0 >= q0)
345 myprow0 = mypcol0 = -1;
346 assert((myprow0 < p0 && mypcol0 < q0) || (myprow0 == -1 && mypcol0 == -1));
347 Cblacs_gridinfo(mb->ctxt, &p1, &q1, &myprow1, &mypcol1);
348 if (myprow1 >= p1 || mypcol1 >= q1)
349 myprow1 = mypcol1 = -1;
350 assert((myprow1 < p1 && mypcol1 < q1) || (myprow1 == -1 && mypcol1 == -1));
351 /* exchange the missing parameters among the processors: shape of grids and
352 * location of the processors */
353 param = (Int *) mr2d_malloc(3 * (nprocs * 2 + NBPARAM) * sizeof(Int));
354 ra = param + nprocs * 2 + NBPARAM;
355 ca = param + (nprocs * 2 + NBPARAM) * 2;
356 for (i = 0; i < nprocs * 2 + NBPARAM; i++)
357 param[i] = MAGIC_MAX;
358 proc0 = param + NBPARAM;
359 proc1 = param + NBPARAM + nprocs;
360 /* we calulate proc0 and proc1 that will give the number of a proc in
361 * respectively a or b in the global context */
362 if (myprow0 >= 0) {
363 proc0[myprow0 * q0 + mypcol0] = mypnum;
364 param[0] = p0;
365 param[1] = q0;
366 param[4] = ma->m;
367 param[5] = ma->n;
368 param[6] = ma->nbrow;
369 param[7] = ma->nbcol;
370 param[8] = ma->sprow;
371 param[9] = ma->spcol;
372 param[10] = ia;
373 param[11] = ja;
374 }
375 if (myprow1 >= 0) {
376 proc1[myprow1 * q1 + mypcol1] = mypnum;
377 param[2] = p1;
378 param[3] = q1;
379 param[12] = mb->m;
380 param[13] = mb->n;
381 param[14] = mb->nbrow;
382 param[15] = mb->nbcol;
383 param[16] = mb->sprow;
384 param[17] = mb->spcol;
385 param[18] = ib;
386 param[19] = jb;
387 }
388 Cigamn2d(gcontext, "All", "H", 2 * nprocs + NBPARAM, (Int)1, param, 2 * nprocs + NBPARAM,
389 ra, ca, 2 * nprocs + NBPARAM, (Int)-1, (Int)-1);
390 newa = *ma;
391 newb = *mb;
392 ma = &newa;
393 mb = &newb;
394 if (myprow0 == -1) {
395 p0 = param[0];
396 q0 = param[1];
397 ma->m = param[4];
398 ma->n = param[5];
399 ma->nbrow = param[6];
400 ma->nbcol = param[7];
401 ma->sprow = param[8];
402 ma->spcol = param[9];
403 ia = param[10];
404 ja = param[11];
405 }
406 if (myprow1 == -1) {
407 p1 = param[2];
408 q1 = param[3];
409 mb->m = param[12];
410 mb->n = param[13];
411 mb->nbrow = param[14];
412 mb->nbcol = param[15];
413 mb->sprow = param[16];
414 mb->spcol = param[17];
415 ib = param[18];
416 jb = param[19];
417 }
418 for (i = 0; i < NBPARAM; i++) {
419 if (param[i] == MAGIC_MAX) {
420 fprintf(stderr, "xxGEMR2D:something wrong in the parameters\n");
421 exit(1);
422 }
423 }
424#ifndef NDEBUG
425 for (i = 0; i < p0 * q0; i++)
426 assert(proc0[i] >= 0 && proc0[i] < nprocs);
427 for (i = 0; i < p1 * q1; i++)
428 assert(proc1[i] >= 0 && proc1[i] < nprocs);
429#endif
430 /* check the validity of the parameters */
431 paramcheck(ma, ia, ja, m, n, p0, q0, gcontext);
432 paramcheck(mb, ib, jb, m, n, p1, q1, gcontext);
433 /* we change the problem so that ia < a->nbrow ... andia + m = a->m ... */
434 {
435 Int decal;
436 ia = changeorigin(myprow0, ma->sprow, p0,
437 ma->nbrow, ia, &decal, &ma->sprow);
438 ptrmyblock += decal;
439 ja = changeorigin(mypcol0, ma->spcol, q0,
440 ma->nbcol, ja, &decal, &ma->spcol);
441 ptrmyblock += decal * ma->lda;
442 ma->m = ia + m;
443 ma->n = ja + n;
444 ib = changeorigin(myprow1, mb->sprow, p1,
445 mb->nbrow, ib, &decal, &mb->sprow);
446 ptrmynewblock += decal;
447 jb = changeorigin(mypcol1, mb->spcol, q1,
448 mb->nbcol, jb, &decal, &mb->spcol);
449 ptrmynewblock += decal * mb->lda;
450 mb->m = ib + m;
451 mb->n = jb + n;
452 if (p0 == 1)
453 ma->nbrow = ma->m;
454 if (q0 == 1)
455 ma->nbcol = ma->n;
456 if (p1 == 1)
457 mb->nbrow = mb->m;
458 if (q1 == 1)
459 mb->nbcol = mb->n;
460#ifndef NDEBUG
461 paramcheck(ma, ia, ja, m, n, p0, q0, gcontext);
462 paramcheck(mb, ib, jb, m, n, p1, q1, gcontext);
463#endif
464 }
465 /* We compute the size of the memory buffer ( we choose the worst case,
466 * when the buffer sizes == the memory block sizes). */
467 if (myprow0 >= 0 && mypcol0 >= 0) {
468 /* Initialize pointer variables */
469 setmemory(&ptrsendbuff, memoryblocksize(ma));
470 }; /* if (mypnum < p0 * q0) */
471 if (myprow1 >= 0 && mypcol1 >= 0) {
472 /* Initialize pointer variables */
473 setmemory(&ptrrecvbuff, memoryblocksize(mb));
474 }; /* if (mypnum < p1 * q1) */
475 /* allocing room for the tabs, alloc for the worst case,local_n or local_m
476 * intervals, in fact the worst case should be less, perhaps half that,I
477 * should think of that one day. */
478 h_inter = (IDESC *) mr2d_malloc(DIVUP(ma->n, q0 * ma->nbcol) *
479 ma->nbcol * sizeof(IDESC));
480 v_inter = (IDESC *) mr2d_malloc(DIVUP(ma->m, p0 * ma->nbrow)
481 * ma->nbrow * sizeof(IDESC));
482 /* We go for the scanning of indices. For each processor including mypnum,
483 * we fill the sendbuff buffer (scanD0(SENDBUFF)) and when it is done send
484 * it. Then for each processor, we compute the size of message to be
485 * receive scanD0(SIZEBUFF)), post a receive and then allocate the elements
486 * of recvbuff the right place (scanD)(RECVBUFF)) */
487 recvptr = ptrrecvbuff;
488 {
489 Int tot, myrang, step, sens;
490 Int *sender, *recver;
491 Int mesending, merecving;
492 tot = max(p0 * q0, p1 * q1);
493 init_chenille(mypnum, nprocs, p0 * q0, proc0, p1 * q1, proc1,
494 &sender, &recver, &myrang);
495 if (myrang == -1)
496 goto after_comm;
497 mesending = myprow0 >= 0;
498 assert(sender[myrang] >= 0 || !mesending);
499 assert(!mesending || proc0[sender[myrang]] == mypnum);
500 merecving = myprow1 >= 0;
501 assert(recver[myrang] >= 0 || !merecving);
502 assert(!merecving || proc1[recver[myrang]] == mypnum);
503 step = tot - 1 - myrang;
504 do {
505 for (sens = 0; sens < 2; sens++) {
506 /* be careful here, when we communicating with ourselves, we must
507 * send first (myrang > step == 0) */
508 if (mesending && recver[step] >= 0 &&
509 (sens == 0)) {
510 i = recver[step] / q1;
511 j = recver[step] % q1;
512 vinter_nb = scan_intervals('r', ia, ib, m, ma, mb, p0, p1, myprow0, i,
513 v_inter);
514 hinter_nb = scan_intervals('c', ja, jb, n, ma, mb, q0, q1, mypcol0, j,
515 h_inter);
516 scanD0(uplo, diag, SENDBUFF, ptrsendbuff, &sendsize,
517 m, n, ma, ia, ja, p0, q0, mb, ib, jb, p1, q1,
518 v_inter, vinter_nb, h_inter, hinter_nb,
519 ptrmyblock);
520 } /* if (mesending...) { */
521 if (mesending && recver[step] >= 0 &&
522 (sens == myrang > step)) {
523 i = recver[step] / q1;
524 j = recver[step] % q1;
525 if (sendsize > 0
526 && (step != myrang || !merecving)
527 ) {
528 Csgesd2d(gcontext, sendsize, (Int)1, ptrsendbuff, sendsize,
529 (Int)0, proc1[i * q1 + j]);
530 } /* sendsize > 0 */
531 } /* if (mesending ... */
532 if (merecving && sender[step] >= 0 &&
533 (sens == myrang <= step)) {
534 i = sender[step] / q0;
535 j = sender[step] % q0;
536 vinter_nb = scan_intervals('r', ib, ia, m, mb, ma, p1, p0, myprow1, i,
537 v_inter);
538 hinter_nb = scan_intervals('c', jb, ja, n, mb, ma, q1, q0, mypcol1, j,
539 h_inter);
540 scanD0(uplo, diag, SIZEBUFF, ptrNULL, &recvsize,
541 m, n, ma, ia, ja, p0, q0, mb, ib, jb, p1, q1,
542 v_inter, vinter_nb, h_inter, hinter_nb, ptrNULL);
543 if (recvsize > 0) {
544 if (step == myrang && mesending) {
545 Clacpy(recvsize, (Int)1,
546 ptrsendbuff, recvsize,
547 ptrrecvbuff, recvsize);
548 } else {
549 Csgerv2d(gcontext, recvsize, (Int)1, ptrrecvbuff, recvsize,
550 (Int)0, proc0[i * q0 + j]);
551 }
552 } /* recvsize > 0 */
553 } /* if (merecving ...) */
554 if (merecving && sender[step] >= 0 && sens == 1) {
555 scanD0(uplo, diag, RECVBUFF, ptrrecvbuff, &recvsize,
556 m, n, ma, ia, ja, p0, q0, mb, ib, jb, p1, q1,
557 v_inter, vinter_nb, h_inter, hinter_nb, ptrmynewblock);
558 } /* if (merecving...) */
559 } /* for (sens = 0) */
560 step -= 1;
561 if (step < 0)
562 step = tot - 1;
563 } while (step != tot - 1 - myrang);
564after_comm:
565 free(sender);
566 } /* { int tot,nr,ns ...} */
567 /* don't forget to clean up things! */
568 if (myprow1 >= 0 && mypcol1 >= 0) {
569 freememory((char *) ptrrecvbuff);
570 };
571 if (myprow0 >= 0 && mypcol0 >= 0) {
572 freememory((char *) ptrsendbuff);
573 };
574 if (nprow != 1)
575 Cblacs_gridexit(gcontext);
576 free(v_inter);
577 free(h_inter);
578 free(param);
579}/* distrib */
581init_chenille(Int mypnum, Int nprocs, Int n0, Int *proc0, Int n1, Int *proc1, Int **psend, Int **precv, Int *myrang)
582{
583 Int ns, nr, i, tot;
584 Int *sender, *recver, *g0, *g1;
585 tot = max(n0, n1);
586 sender = (Int *) mr2d_malloc((nprocs + tot) * sizeof(Int) * 2);
587 recver = sender + tot;
588 *psend = sender;
589 *precv = recver;
590 g0 = recver + tot;
591 g1 = g0 + nprocs;
592 for (i = 0; i < nprocs; i++) {
593 g0[i] = -1;
594 g1[i] = -1;
595 }
596 for (i = 0; i < tot; i++) {
597 sender[i] = -1;
598 recver[i] = -1;
599 }
600 for (i = 0; i < n0; i++)
601 g0[proc0[i]] = i;
602 for (i = 0; i < n1; i++)
603 g1[proc1[i]] = i;
604 ns = 0;
605 nr = 0;
606 *myrang = -1;
607 for (i = 0; i < nprocs; i++)
608 if (g0[i] >= 0 && g1[i] >= 0) {
609 if (i == mypnum)
610 *myrang = nr;
611 sender[ns] = g0[i];
612 ns += 1;
613 recver[nr] = g1[i];
614 nr += 1;
615 assert(ns <= n0 && nr <= n1 && nr == ns);
616 }
617 for (i = 0; i < nprocs; i++)
618 if (g0[i] >= 0 && g1[i] < 0) {
619 if (i == mypnum)
620 *myrang = ns;
621 sender[ns] = g0[i];
622 ns += 1;
623 assert(ns <= n0);
624 }
625 for (i = 0; i < nprocs; i++)
626 if (g1[i] >= 0 && g0[i] < 0) {
627 if (i == mypnum)
628 *myrang = nr;
629 recver[nr] = g1[i];
630 nr += 1;
631 assert(nr <= n1);
632 }
633}
634void
635Clacpy(Int m, Int n, float *a, Int lda, float *b, Int ldb)
636{
637 Int i, j;
638 lda -= m;
639 ldb -= m;
640 assert(lda >= 0 && ldb >= 0);
641 for (j = 0; j < n; j++) {
642 for (i = 0; i < m; i++)
643 *b++ = *a++;
644 b += ldb;
645 a += lda;
646 }
647}
649gridreshape(Int *ctxtp)
650{
651 Int ori, final; /* original context, and new context created, with
652 * line form */
653 Int nprow, npcol, myrow, mycol;
654 Int *usermap;
655 Int i, j;
656 ori = *ctxtp;
657 Cblacs_gridinfo(ori, &nprow, &npcol, &myrow, &mycol);
658 usermap = mr2d_malloc(sizeof(Int) * nprow * npcol);
659 for (i = 0; i < nprow; i++)
660 for (j = 0; j < npcol; j++) {
661 usermap[i + j * nprow] = Cblacs_pnum(ori, i, j);
662 }
663 /* Cblacs_get(0, 0, &final); */
664 Cblacs_get(ori, (Int)10, &final);
665 Cblacs_gridmap(&final, usermap, (Int)1, (Int)1, nprow * npcol);
666 *ctxtp = final;
667 free(usermap);
668}
#define Int
Definition Bconfig.h:22
void Cblacs_gridmap()
static2 Int inter_len()
Int memoryblocksize()
Int changeorigin()
#define freememory
Definition pstrmr.c:235
#define scan_intervals
Definition pstrmr.c:236
void Csgerv2d()
void Cpstrmr2do()
void Cblacs_gridexit()
#define max(A, B)
Definition pstrmr.c:192
void Cigsum2d()
#define static2
Definition pstrmr.c:158
#define scanD0
Definition pstrmr.c:232
Int Cblacs_pnum()
static2 void gridreshape()
#define DIVUP(a, b)
Definition pstrmr.c:194
#define SIZEBUFF
Definition pstrmr.c:247
#define fortran_mr2dnew
Definition pstrmr.c:169
#define Clacpy
Definition pstrmr.c:173
void Cigamx2d()
void Cblacs_pinfo()
#define NBPARAM
Definition pstrmr.c:301
void Cblacs_pcoord()
void Cigamn2d()
Int localsize()
#define setmemory
Definition pstrmr.c:234
void Cblacs_get()
void Cigerv2d()
void Cigebs2d()
#define dispmat
Definition pstrmr.c:233
static2 void buff2block()
#define SENDBUFF
Definition pstrmr.c:245
void paramcheck()
void Cblacs_setup()
#define RECVBUFF
Definition pstrmr.c:246
static2 Int block2buff()
void Cblacs_gridinit()
void Cblacs_gridinfo()
#define MAGIC_MAX
Definition pstrmr.c:302
void Csetpvmtids()
Int localindice()
void Cigesd2d()
#define slacpy_
Definition pstrmr.c:171
void Cpstrmr2d()
#define DESCLEN
Definition pstrmr.c:257
Int ppcm()
void Cigebr2d()
void Cblacs_exit()
#define fortran_mr2d
Definition pstrmr.c:168
void * mr2d_malloc()
void Csgesd2d()
static2 void init_chenille()
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 n
Definition pcgemr.c:167
Int lda
Definition pcgemr.c:172