1 #include "redist.h"
2 /* $Id: psgemr2.c,v 1.1.1.1 2000/02/15 18:04:09 susan Exp $
3 *
4 * some functions used by the psgemr2d routine see file psgemr.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 psgemr2do_
11 #define fortran_mr2dnew psgemr2d_
12 #elif defined(UpCase)
13 #define fortran_mr2dnew PSGEMR2D
14 #define fortran_mr2d PSGEMR2DO
15 #define scopy_ SCOPY
16 #define slacpy_ SLACPY
17 #else
18 #define fortran_mr2d psgemr2do
19 #define fortran_mr2dnew psgemr2d
20 #define scopy_ scopy
21 #define slacpy_ slacpy
22 #endif
23 #define Clacpy Csgelacpy
24 void Clacpy();
25 typedef struct {
26 int desctype;
27 int ctxt;
28 int m;
29 int n;
30 int nbrow;
31 int nbcol;
32 int sprow;
33 int spcol;
34 int lda;
35 } MDESC;
36 #define BLOCK_CYCLIC_2D 1
37 typedef struct {
38 int lstart;
39 int len;
40 } IDESC;
41 #define SHIFT(row,sprow,nbrow) ((row)-(sprow)+ ((row) >= (sprow) ? 0 : (nbrow)))
42 #define max(A,B) ((A)>(B)?(A):(B))
43 #define min(A,B) ((A)>(B)?(B):(A))
44 #define DIVUP(a,b) ( ((a)-1) /(b)+1)
45 #define ROUNDUP(a,b) (DIVUP(a,b)*(b))
46 #ifdef MALLOCDEBUG
47 #define malloc mymalloc
48 #define free myfree
49 #define realloc myrealloc
50 #endif
51 /* Cblacs */
52 extern void Cblacs_pcoord();
53 extern int Cblacs_pnum();
54 extern void Csetpvmtids();
55 extern void Cblacs_get();
56 extern void Cblacs_pinfo();
57 extern void Cblacs_gridinfo();
58 extern void Cblacs_gridinit();
59 extern void Cblacs_exit();
60 extern void Cblacs_gridexit();
61 extern void Cblacs_setup();
62 extern void Cigebs2d();
63 extern void Cigebr2d();
64 extern void Cigesd2d();
65 extern void Cigerv2d();
66 extern void Cigsum2d();
67 extern void Cigamn2d();
68 extern void Cigamx2d();
69 extern void Csgesd2d();
70 extern void Csgerv2d();
71 /* lapack */
72 void slacpy_();
73 /* aux fonctions */
74 extern int localindice();
75 extern void *mr2d_malloc();
76 extern int ppcm();
77 extern int localsize();
78 extern int memoryblocksize();
79 extern int changeorigin();
80 extern void paramcheck();
81 /* tools and others function */
82 #define scanD0 sgescanD0
83 #define dispmat sgedispmat
84 #define setmemory sgesetmemory
85 #define freememory sgefreememory
86 #define scan_intervals sgescan_intervals
87 extern void scanD0();
88 extern void dispmat();
89 extern void setmemory();
90 extern void freememory();
91 extern int scan_intervals();
92 extern void Cpsgemr2do();
93 extern void Cpsgemr2d();
94 /* some defines for Cpsgemr2do */
95 #define SENDBUFF 0
96 #define RECVBUFF 1
97 #define SIZEBUFF 2
98 #if 0
99 #define DEBUG
100 #endif
101 #ifndef DEBUG
102 #define NDEBUG
103 #endif
104 #include <stdio.h>
105 #include <stdlib.h>
106 #include <string.h>
107 #include <assert.h>
108 #include <ctype.h>
109 /* Created March 1993 by B. Tourancheau (See sccs for modifications). */
110 /************************************************************************/
111 /* Set the memory space with the malloc function */
112 void
setmemory(adpointer,blocksize)113 setmemory(adpointer, blocksize)
114 float **adpointer;
115 int blocksize;
116 {
117 assert(blocksize >= 0);
118 if (blocksize == 0) {
119 *adpointer = NULL;
120 return;
121 }
122 *adpointer = (float *) mr2d_malloc(
123 blocksize * sizeof(float));
124 }
125 /******************************************************************/
126 /* Free the memory space after the malloc */
127 void
freememory(ptrtobefreed)128 freememory(ptrtobefreed)
129 float *ptrtobefreed;
130 {
131 if (ptrtobefreed == NULL)
132 return;
133 free((char *) ptrtobefreed);
134 }
135 /* extern functions for intersect() extern scopy_(); */
136 /* scan_intervals: scans two distributions in one dimension, and compute the
137 * intersections on the local processor. result must be long enough to
138 * contains the result that are stocked in IDESC structure, the function
139 * returns the number of intersections found */
140 int
scan_intervals(type,ja,jb,n,ma,mb,q0,q1,col0,col1,result)141 scan_intervals(type, ja, jb, n, ma, mb, q0, q1, col0, col1,
142 result)
143 char type;
144 int ja, jb, n, q0, q1, col0, col1;
145 MDESC *ma, *mb;
146 IDESC *result;
147 {
148 int offset, j0, j1, templatewidth0, templatewidth1, nbcol0, nbcol1;
149 int l; /* local indice on the beginning of the interval */
150 assert(type == 'c' || type == 'r');
151 nbcol0 = (type == 'c' ? ma->nbcol : ma->nbrow);
152 nbcol1 = (type == 'c' ? mb->nbcol : mb->nbrow);
153 templatewidth0 = q0 * nbcol0;
154 templatewidth1 = q1 * nbcol1;
155 {
156 int sp0 = (type == 'c' ? ma->spcol : ma->sprow);
157 int sp1 = (type == 'c' ? mb->spcol : mb->sprow);
158 j0 = SHIFT(col0, sp0, q0) * nbcol0 - ja;
159 j1 = SHIFT(col1, sp1, q1) * nbcol1 - jb;
160 }
161 offset = 0;
162 l = 0;
163 /* a small check to verify that the submatrix begin inside the first block
164 * of the original matrix, this done by a sort of coordinate change at the
165 * beginning of the Cpsgemr2d */
166 assert(j0 + nbcol0 > 0);
167 assert(j1 + nbcol1 > 0);
168 while ((j0 < n) && (j1 < n)) {
169 int end0, end1;
170 int start, end;
171 end0 = j0 + nbcol0;
172 end1 = j1 + nbcol1;
173 if (end0 <= j1) {
174 j0 += templatewidth0;
175 l += nbcol0;
176 continue;
177 }
178 if (end1 <= j0) {
179 j1 += templatewidth1;
180 continue;
181 }
182 /* compute the raw intersection */
183 start = max(j0, j1);
184 start = max(start, 0);
185 /* the start is correct now, update the corresponding fields */
186 result[offset].lstart = l + start - j0;
187 end = min(end0, end1);
188 if (end0 == end) {
189 j0 += templatewidth0;
190 l += nbcol0;
191 }
192 if (end1 == end)
193 j1 += templatewidth1;
194 /* throw the limit if they go out of the matrix */
195 end = min(end, n);
196 assert(end > start);
197 /* it is a bit tricky to see why the length is always positive after all
198 * this min and max, first we have the property that every interval
199 * considered is at least partly into the submatrix, second we arrive
200 * here only if the raw intersection is non-void, if we remove a limit
201 * that means the corresponding frontier is in both intervals which
202 * proove the final interval is non-void, clear ?? */
203 result[offset].len = end - start;
204 offset += 1;
205 } /* while */
206 return offset;
207 }
208