ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pctrmr2.c
Go to the documentation of this file.
1 #include "redist.h"
2 /* $Id: pctrmr2.c,v 1.1.1.1 2000/02/15 18:04:10 susan Exp $
3  *
4  * some functions used by the pctrmr2d routine see file pctrmr.c for more
5  * documentation.
6  *
7  * Created March 1993 by B. Tourancheau (See sccs for modifications). */
8 #define static2 static
9 #if defined(Add_) || defined(f77IsF2C)
10 #define fortran_mr2d pctrmr2do_
11 #define fortran_mr2dnew pctrmr2d_
12 #elif defined(UpCase)
13 #define fortran_mr2dnew PCTRMR2D
14 #define fortran_mr2d PCTRMR2DO
15 #define ccopy_ CCOPY
16 #define clacpy_ CLACPY
17 #else
18 #define fortran_mr2d pctrmr2do
19 #define fortran_mr2dnew pctrmr2d
20 #define ccopy_ ccopy
21 #define clacpy_ clacpy
22 #endif
23 #define Clacpy Cctrlacpy
24 void Clacpy();
25 typedef struct {
26  float r, i;
27 } complex;
28 typedef struct {
29  int desctype;
30  int ctxt;
31  int m;
32  int n;
33  int nbrow;
34  int nbcol;
35  int sprow;
36  int spcol;
37  int lda;
38 } MDESC;
39 #define BLOCK_CYCLIC_2D 1
40 typedef struct {
41  int gstart;
42  int len;
43 } IDESC;
44 #define SHIFT(row,sprow,nbrow) ((row)-(sprow)+ ((row) >= (sprow) ? 0 : (nbrow)))
45 #define max(A,B) ((A)>(B)?(A):(B))
46 #define min(A,B) ((A)>(B)?(B):(A))
47 #define DIVUP(a,b) ( ((a)-1) /(b)+1)
48 #define ROUNDUP(a,b) (DIVUP(a,b)*(b))
49 #ifdef MALLOCDEBUG
50 #define malloc mymalloc
51 #define free myfree
52 #define realloc myrealloc
53 #endif
54 /* Cblacs */
55 extern void Cblacs_pcoord();
56 extern int Cblacs_pnum();
57 extern void Csetpvmtids();
58 extern void Cblacs_get();
59 extern void Cblacs_pinfo();
60 extern void Cblacs_gridinfo();
61 extern void Cblacs_gridinit();
62 extern void Cblacs_exit();
63 extern void Cblacs_gridexit();
64 extern void Cblacs_setup();
65 extern void Cigebs2d();
66 extern void Cigebr2d();
67 extern void Cigesd2d();
68 extern void Cigerv2d();
69 extern void Cigsum2d();
70 extern void Cigamn2d();
71 extern void Cigamx2d();
72 extern void Ccgesd2d();
73 extern void Ccgerv2d();
74 /* lapack */
75 void clacpy_();
76 /* aux fonctions */
77 extern int localindice();
78 extern void *mr2d_malloc();
79 extern int ppcm();
80 extern int localsize();
81 extern int memoryblocksize();
82 extern int changeorigin();
83 extern void paramcheck();
84 /* tools and others function */
85 #define scanD0 ctrscanD0
86 #define dispmat ctrdispmat
87 #define setmemory ctrsetmemory
88 #define freememory ctrfreememory
89 #define scan_intervals ctrscan_intervals
90 extern void scanD0();
91 extern void dispmat();
92 extern void setmemory();
93 extern void freememory();
94 extern int scan_intervals();
95 extern void Cpctrmr2do();
96 extern void Cpctrmr2d();
97 /* some defines for Cpctrmr2do */
98 #define SENDBUFF 0
99 #define RECVBUFF 1
100 #define SIZEBUFF 2
101 #if 0
102 #define DEBUG
103 #endif
104 #ifndef DEBUG
105 #define NDEBUG
106 #endif
107 #include <stdio.h>
108 #include <stdlib.h>
109 #include <string.h>
110 #include <assert.h>
111 #include <ctype.h>
112 /* Created March 1993 by B. Tourancheau (See sccs for modifications). */
113 /************************************************************************/
114 /* Set the memory space with the malloc function */
115 void
116 setmemory(adpointer, blocksize)
117  complex **adpointer;
118  int blocksize;
119 {
120  assert(blocksize >= 0);
121  if (blocksize == 0) {
122  *adpointer = NULL;
123  return;
124  }
125  *adpointer = (complex *) mr2d_malloc(
126  blocksize * sizeof(complex));
127 }
128 /******************************************************************/
129 /* Free the memory space after the malloc */
130 void
131 freememory(ptrtobefreed)
132  complex *ptrtobefreed;
133 {
134  if (ptrtobefreed == NULL)
135  return;
136  free((char *) ptrtobefreed);
137 }
138 /* extern functions for intersect() extern ccopy_(); */
139 /**************************************************************/
140 /* return the number of elements int the column after i and the distance of
141  * the first one from i, i,j can be negative out of borns, the number of
142  * elements returned can be negative (means 0) */
143 static2 int
144 insidemat(uplo, diag, i, j, m, n, offset)
145  int m, n, i, j; /* coordonnees de depart, taille de la sous-matrice */
146  char *uplo, *diag;
147  int *offset;
148 {
149  /* tests outside mxn */
150  assert(j >= 0 && j < n);
151  assert(i >= 0);
152  if (toupper(*uplo) == 'U') {
153  int nbline; /* number of lines in the j_th column */
154  int virtualnbline; /* number of line if we were not limited by m */
155  *offset = 0;
156  virtualnbline = max(m - n, 0) + j + (toupper(*diag) == 'N');
157  nbline = min(virtualnbline, m);
158  return nbline - i;
159  } else {
160  int firstline; /* first line in the j_th column */
161  int diagcol; /* column where the diag begin */
162  int virtualline; /* virtual first line if the matrix was extended with
163  * negative indices */
164  int off;
165  diagcol = max(n - m, 0);;
166  virtualline = j - diagcol + (toupper(*diag) == 'U');
167  firstline = max(0, virtualline);
168  off = max(firstline - i, 0);
169  *offset = off;
170  i += off;
171  return m - i;
172  }
173 }/* insidemat() */
174 /********************************************************************/
175 /* Execute an action on the local memories when an intersection occurs (the
176  * action can be the filling of the memory buffer, the count of the memory
177  * buffer size or the setting of the memory with the element received) */
178 static2 void
179 intersect(uplo, diag,
180  j, start, end,
181  action,
182  ptrsizebuff, pptrbuff, ptrblock,
183  m, n,
184  ma, ia, ja, templateheight0, templatewidth0,
185  mb, ib, jb, templateheight1, templatewidth1)
186  int action, *ptrsizebuff;
187  int j, start, end;
188  complex **pptrbuff, *ptrblock;
189  int templateheight0, templatewidth0;
190  int templateheight1, templatewidth1;
191  MDESC *ma, *mb;
192  int ia, ja, ib, jb, m, n;
193  char *uplo, *diag;
194 /* Execute the action on the local memory for the current interval and
195  * increment pptrbuff and ptrsizebuff of the intervalsize */
196 /* Notice that if the interval is contigous in the virtual matrice, it is
197  * also contigous in the real one ! */
198 {
199  /* int un = 1; only when we use dcopy instead of memcpy */
200  complex *ptrstart;
201  int offset, nbline;
202  int intervalsize;
203  assert(start < end);
204  assert(j >= 0 && j < n);
205  nbline =
206  insidemat(uplo, diag, start, j, m, n, &offset);
207  if (nbline <= 0)
208  return;
209  start += offset;
210  if (start >= end)
211  return;
212  intervalsize = min(end - start, nbline);
213  (*ptrsizebuff) += intervalsize;
214  switch (action) {
215  case SENDBUFF: /* fill buff with local elements to be sent */
216  ptrstart = ptrblock + localindice(start + ia, j + ja,
217  templateheight0, templatewidth0, ma);
218  memcpy((char *) (*pptrbuff), (char *) ptrstart,
219  intervalsize * sizeof(complex));
220  /* ccopy_(&intervalsize, (char *) (ptrstart), &un, (char *) (*pptrbuff),
221  * &un); */
222  (*pptrbuff) += intervalsize;
223  break;
224  case RECVBUFF: /* fill local memory with the values received */
225  ptrstart = ptrblock + localindice(start + ib, j + jb,
226  templateheight1, templatewidth1, mb);
227  memcpy((char *) ptrstart, (char *) (*pptrbuff),
228  intervalsize * sizeof(complex));
229  /* ccopy_(&intervalsize, (char *) (*pptrbuff), &un, (char *) (ptrstart),
230  * &un); */
231  (*pptrbuff) += intervalsize;
232  break;
233  case SIZEBUFF: /* computation of sizebuff */
234  break;
235  default:
236  printf("action is %d outside the scope of the case [0..2] !! \n ", action);
237  exit(0);
238  break;
239  }; /* switch (action) */
240 }/* intersect() */
241 /* scan_intervals: scans two distributions in one dimension, and compute the
242  * intersections on the local processor. result must be long enough to
243  * contains the result that are stocked in IDESC structure, the function
244  * returns the number of intersections found */
245 int
246 scan_intervals(type, ja, jb, n, ma, mb, q0, q1, col0, col1,
247  result)
248  char type;
249  int ja, jb, n, q0, q1, col0, col1;
250  MDESC *ma, *mb;
251  IDESC *result;
252 {
253  int offset, j0, j1, templatewidth0, templatewidth1, nbcol0, nbcol1;
254  int l; /* local indice on the beginning of the interval */
255  assert(type == 'c' || type == 'r');
256  nbcol0 = (type == 'c' ? ma->nbcol : ma->nbrow);
257  nbcol1 = (type == 'c' ? mb->nbcol : mb->nbrow);
258  templatewidth0 = q0 * nbcol0;
259  templatewidth1 = q1 * nbcol1;
260  {
261  int sp0 = (type == 'c' ? ma->spcol : ma->sprow);
262  int sp1 = (type == 'c' ? mb->spcol : mb->sprow);
263  j0 = SHIFT(col0, sp0, q0) * nbcol0 - ja;
264  j1 = SHIFT(col1, sp1, q1) * nbcol1 - jb;
265  }
266  offset = 0;
267  l = 0;
268  /* a small check to verify that the submatrix begin inside the first block
269  * of the original matrix, this done by a sort of coordinate change at the
270  * beginning of the Cpctrmr2d */
271  assert(j0 + nbcol0 > 0);
272  assert(j1 + nbcol1 > 0);
273  while ((j0 < n) && (j1 < n)) {
274  int end0, end1;
275  int start, end;
276  end0 = j0 + nbcol0;
277  end1 = j1 + nbcol1;
278  if (end0 <= j1) {
279  j0 += templatewidth0;
280  l += nbcol0;
281  continue;
282  }
283  if (end1 <= j0) {
284  j1 += templatewidth1;
285  continue;
286  }
287  /* compute the raw intersection */
288  start = max(j0, j1);
289  start = max(start, 0);
290  /* the start is correct now, update the corresponding fields */
291  result[offset].gstart = start;
292  end = min(end0, end1);
293  if (end0 == end) {
294  j0 += templatewidth0;
295  l += nbcol0;
296  }
297  if (end1 == end)
298  j1 += templatewidth1;
299  /* throw the limit if they go out of the matrix */
300  end = min(end, n);
301  assert(end > start);
302  /* it is a bit tricky to see why the length is always positive after all
303  * this min and max, first we have the property that every interval
304  * considered is at least partly into the submatrix, second we arrive
305  * here only if the raw intersection is non-void, if we remove a limit
306  * that means the corresponding frontier is in both intervals which
307  * proove the final interval is non-void, clear ?? */
308  result[offset].len = end - start;
309  offset += 1;
310  } /* while */
311  return offset;
312 }
313 /*********************************************************************/
314 /* Do the scanning of intervals and the requested action */
315 void
316 scanD0(uplo, diag, action, ptrbuff, ptrsizebuff,
317  m, n,
318  ma, ia, ja, p0, q0,
319  mb, ib, jb, p1, q1,
320  v_inter, vinter_nb,
321  h_inter, hinter_nb,
322  ptrblock)
323  int action, /* # of the action done on the intersected intervals */
324  *ptrsizebuff; /* size of the communication ptrbuffer (chosen to be
325  * an output parameter in every cases) */
326  complex *ptrbuff /* address of the communication ptrbuffer (a
327  * suffisant memory space is supposed to be allocated
328  before the call) */ , *ptrblock;
329  int p0, q0, p1, q1;
330  IDESC *v_inter, *h_inter;
331  int vinter_nb, hinter_nb;
332  int m, n;
333  int ia, ja, ib, jb;
334  MDESC *ma, *mb;
335  char *uplo, *diag;
336 {/* Rmk: the a+au type addresses are strict bounds as a+au does not belong to
337  * the [a..a+au-1] interval of length au */
338  int templateheight1, templatewidth1;
339  int templateheight0, templatewidth0;
340  int h, v; /* for scanning the intervals */
341  /* initializations */
342  templateheight1 = p1 * mb->nbrow;
343  templateheight0 = p0 * ma->nbrow;
344  templatewidth1 = q1 * mb->nbcol;
345  templatewidth0 = q0 * ma->nbcol;
346  /* we now will deal will logical grids, that's to say we change our
347  * numbering of processors so that (0,0) begin on logical processor (0,0) */
348  /* in case we will not enter the while loop */
349  (*ptrsizebuff) = 0;
350  for (h = 0; h < hinter_nb; h++)
351  for (v = 0; v < vinter_nb; v++) {
352  int j;
353  for (j = 0; j < h_inter[h].len; j++)
354  intersect(uplo, diag, j + h_inter[h].gstart,
355  v_inter[v].gstart, v_inter[v].gstart + v_inter[v].len,
356  action, ptrsizebuff, &ptrbuff, ptrblock, m, n,
357  ma, ia, ja, templateheight0, templatewidth0,
358  mb, ib, jb, templateheight1, templatewidth1);
359  }
360 }/* scanD0() */
intersect
static2 void intersect(char *uplo, char *diag, int j, int start, int end, int action, int *ptrsizebuff, complex **pptrbuff, complex *ptrblock, int m, int n, MDESC *ma, int ia, int ja, int templateheight0, int templatewidth0, MDESC *mb, int ib, int jb, int templateheight1, int templatewidth1)
Definition: pctrmr2.c:179
setmemory
#define setmemory
Definition: pctrmr2.c:87
dispmat
#define dispmat
Definition: pctrmr2.c:86
SHIFT
#define SHIFT(row, sprow, nbrow)
Definition: pctrmr2.c:44
Cigamx2d
void Cigamx2d()
Cblacs_gridexit
void Cblacs_gridexit()
Cblacs_pnum
int Cblacs_pnum()
Cpctrmr2do
void Cpctrmr2do()
Cigerv2d
void Cigerv2d()
SIZEBUFF
#define SIZEBUFF
Definition: pctrmr2.c:100
MDESC::sprow
int sprow
Definition: pcgemr.c:170
scanD0
#define scanD0
Definition: pctrmr2.c:85
Cblacs_pcoord
void Cblacs_pcoord()
MDESC::nbcol
int nbcol
Definition: pcgemr.c:169
MDESC
Definition: pcgemr.c:163
Cigsum2d
void Cigsum2d()
ppcm
int ppcm()
Cigebs2d
void Cigebs2d()
Cblacs_get
void Cblacs_get()
RECVBUFF
#define RECVBUFF
Definition: pctrmr2.c:99
Cigamn2d
void Cigamn2d()
Cblacs_gridinit
void Cblacs_gridinit()
min
#define min(A, B)
Definition: pctrmr2.c:46
localsize
int localsize()
IDESC
Definition: pcgemr.c:175
IDESC::gstart
int gstart
Definition: pctrmr.c:191
scan_intervals
#define scan_intervals
Definition: pctrmr2.c:89
Cblacs_setup
void Cblacs_setup()
Cpctrmr2d
void Cpctrmr2d()
MDESC::spcol
int spcol
Definition: pcgemr.c:171
Ccgesd2d
void Ccgesd2d()
redist.h
SENDBUFF
#define SENDBUFF
Definition: pctrmr2.c:98
Cblacs_exit
void Cblacs_exit()
Csetpvmtids
void Csetpvmtids()
paramcheck
void paramcheck()
Cblacs_gridinfo
void Cblacs_gridinfo()
freememory
#define freememory
Definition: pctrmr2.c:88
clacpy_
#define clacpy_
Definition: pctrmr2.c:21
localindice
int localindice()
IDESC::len
int len
Definition: pcgemr.c:177
Clacpy
#define Clacpy
Definition: pctrmr2.c:23
Ccgerv2d
void Ccgerv2d()
Cigesd2d
void Cigesd2d()
static2
#define static2
Definition: pctrmr2.c:8
Cblacs_pinfo
void Cblacs_pinfo()
changeorigin
int changeorigin()
Cigebr2d
void Cigebr2d()
max
#define max(A, B)
Definition: pctrmr2.c:45
MDESC::nbrow
int nbrow
Definition: pcgemr.c:168
insidemat
static2 int insidemat(char *uplo, char *diag, int i, int j, int m, int n, int *offset)
Definition: pctrmr2.c:144
mr2d_malloc
void * mr2d_malloc()
complex
Definition: pcgemr.c:160
memoryblocksize
int memoryblocksize()