ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
PB_Cplasca2.c
Go to the documentation of this file.
1 /* ---------------------------------------------------------------------
2 *
3 * -- PBLAS auxiliary routine (version 2.0) --
4 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
5 * and University of California, Berkeley.
6 * April 1, 1998
7 *
8 * ---------------------------------------------------------------------
9 */
10 /*
11 * Include files
12 */
13 #include "../pblas.h"
14 #include "../PBpblas.h"
15 #include "../PBtools.h"
16 #include "../PBblacs.h"
17 #include "../PBblas.h"
18 
19 #ifdef __STDC__
20 void PB_Cplasca2( PBTYP_T * TYPE, char * UPLO, char * CONJUG, int M,
21  int N, char * ALPHA, char * A, int IA, int JA,
22  int * DESCA )
23 #else
24 void PB_Cplasca2( TYPE, UPLO, CONJUG, M, N, ALPHA, A, IA, JA, DESCA )
25 /*
26 * .. Scalar Arguments ..
27 */
28  char * CONJUG, * UPLO;
29  int IA, JA, M, N;
30  char * ALPHA;
31  PBTYP_T * TYPE;
32 /*
33 * .. Array Arguments ..
34 */
35  int * DESCA;
36  char * A;
37 #endif
38 {
39 /*
40 * .. Local Scalars ..
41 */
42  char UploA, herm;
43  int Acol, Arow, Aii, iimax, ilow, imbloc, Aimb1, inbloc, Ainb1,
44  Aoffi, GoEast, GoSouth, ioffd, iupp, izero=0, Ajj, jjmax,
45  Aoffj, joffd, lcmt, lcmt00, Ald, lmbloc, lnbloc, low, lower,
46  m1, Amb, mbloc, mblkd, mblks, Amp, Arcol, Arrow, mycol, myrow,
47  n1, Anb, nbloc, nblkd, nblks, npcol, nprow, Anq, pmb, qnb,
48  size, tmp1, upp, upper;
49  TZSCAL_T scal;
50 /* ..
51 * .. Executable Statements ..
52 *
53 */
54 /*
55 * Quick return if possible
56 */
57  if( ( M <= 0 ) || ( N <= 0 ) ) return;
58 /*
59 * Retrieve process grid information
60 */
61  Cblacs_gridinfo( DESCA[CTXT_], &nprow, &npcol, &myrow, &mycol );
62 /*
63 * Retrieve sub( A )'s local information: Aii, Ajj, Arow, Acol ...
64 */
65  PB_Cainfog2l( M, N, IA, JA, DESCA, nprow, npcol, myrow, mycol, &Aimb1,
66  &Ainb1, &Amp, &Anq, &Aii, &Ajj, &Arow, &Acol, &Arrow, &Arcol );
67 /*
68 * Quick return if I don't own any of sub( A ).
69 */
70  if( ( Amp <= 0 ) || ( Anq <= 0 ) ) return;
71 /*
72 * Initialize lcmt00, mblks, nblks, imbloc, inbloc, lmbloc, lnbloc, ilow, low,
73 * iupp, and upp.
74 */
75  Amb = DESCA[MB_ ]; Anb = DESCA[NB_ ]; Ald = DESCA[LLD_];
76  PB_Cbinfo( 0, Amp, Anq, Aimb1, Ainb1, Amb, Anb, Arrow, Arcol, &lcmt00,
77  &mblks, &nblks, &imbloc, &inbloc, &lmbloc, &lnbloc, &ilow, &low,
78  &iupp, &upp );
79  iimax = ( Aoffi = Aii - 1 ) + ( m1 = Amp );
80  jjmax = ( Aoffj = Ajj - 1 ) + ( n1 = Anq );
81  pmb = ( ( ( Arow < 0 ) || ( nprow == 1 ) ) ? Amb : nprow * Amb );
82  qnb = ( ( ( Acol < 0 ) || ( npcol == 1 ) ) ? Anb : npcol * Anb );
83 
84  UploA = Mupcase( UPLO[0] );
85  upper = ( UploA != CLOWER );
86  lower = ( UploA != CUPPER );
87  herm = ( UploA == CALL ? CNOCONJG : Mupcase( CONJUG[0] ) );
88 
89  size = TYPE->size;
90  scal = ( herm == CCONJG ? TYPE->Fhescal : TYPE->Ftzscal );
91 /*
92 * Handle separately the first row and/or column of the LCM table. Update the
93 * LCM value of the curent block lcmt00, as well as the number of rows and
94 * columns mblks and nblks remaining in the LCM table.
95 */
96  GoSouth = ( lcmt00 > iupp );
97  GoEast = ( lcmt00 < ilow );
98 /*
99 * Go through the table looking for blocks owning diagonal entries.
100 */
101  if( ( !( GoSouth ) ) && ( !( GoEast ) ) )
102  {
103 /*
104 * The upper left block owns diagonal entries lcmt00 >= ilow && lcmt00 <= iupp
105 */
106  scal( C2F_CHAR( UPLO ), &imbloc, &inbloc, &lcmt00, ALPHA,
107  Mptr( A, Aii, Ajj, Ald, size ), &Ald );
108 /*
109 * Decide whether one should go south or east in the table: Go east if
110 * the block below the current one only owns lower entries. If this block,
111 * however, owns diagonals, then go south.
112 */
113  GoSouth = !( GoEast = ( ( lcmt00 - ( iupp - upp + pmb ) ) < ilow ) );
114 
115  if( GoSouth )
116  {
117 /*
118 * When the upper triangular part of sub( A ) should be scaled and one is
119 * planning to go south in the table, it is neccessary to take care of the
120 * remaining columns of these imbloc rows immediately.
121 */
122  if( upper && ( Anq > inbloc ) )
123  {
124  tmp1 = Anq - inbloc;
125  scal( C2F_CHAR( ALL ), &imbloc, &tmp1, &izero, ALPHA,
126  Mptr( A, Aii, Ajj+inbloc, Ald, size ), &Ald );
127  }
128  Aii += imbloc;
129  m1 -= imbloc;
130  }
131  else
132  {
133 /*
134 * When the lower triangular part of sub( A ) should be scaled and one is
135 * planning to go east in the table, it is neccessary to take care of the
136 * remaining rows of these inbloc columns immediately.
137 */
138  if( lower && ( Amp > imbloc ) )
139  {
140  tmp1 = Amp - imbloc;
141  scal( C2F_CHAR( ALL ), &tmp1, &inbloc, &izero, ALPHA,
142  Mptr( A, Aii+imbloc, Ajj, Ald, size ), &Ald );
143  }
144  Ajj += inbloc;
145  n1 -= inbloc;
146  }
147  }
148 
149  if( GoSouth )
150  {
151 /*
152 * Go one step south in the LCM table. Adjust the current LCM value as well as
153 * the local row index in A.
154 */
155  lcmt00 -= ( iupp - upp + pmb ); mblks--; Aoffi += imbloc;
156 /*
157 * While there are blocks remaining that own upper entries, keep going south.
158 * Adjust the current LCM value as well as the local row index in A.
159 */
160  while( ( mblks > 0 ) && ( lcmt00 > upp ) )
161  { lcmt00 -= pmb; mblks--; Aoffi += Amb; }
162 /*
163 * Scale the upper triangular part of sub( A ) we just skipped when necessary.
164 */
165  tmp1 = MIN( Aoffi, iimax ) - Aii + 1;
166  if( upper && ( tmp1 > 0 ) )
167  {
168  scal( C2F_CHAR( ALL ), &tmp1, &n1, &izero, ALPHA,
169  Mptr( A, Aii, Aoffj+1, Ald, size ), &Ald );
170  Aii += tmp1;
171  m1 -= tmp1;
172  }
173 /*
174 * Return if no more row in the LCM table.
175 */
176  if( mblks <= 0 ) return;
177 /*
178 * lcmt00 <= upp. The current block owns either diagonals or lower entries.
179 * Save the current position in the LCM table. After this column has been
180 * completely taken care of, re-start from this row and the next column of
181 * the LCM table.
182 */
183  lcmt = lcmt00; mblkd = mblks; ioffd = Aoffi;
184 
185  mbloc = Amb;
186  while( ( mblkd > 0 ) && ( lcmt >= ilow ) )
187  {
188 /*
189 * A block owning diagonals lcmt00 >= ilow && lcmt00 <= upp has been found.
190 */
191  if( mblkd == 1 ) mbloc = lmbloc;
192  scal( C2F_CHAR( UPLO ), &mbloc, &inbloc, &lcmt, ALPHA,
193  Mptr( A, ioffd+1, Aoffj+1, Ald, size ), &Ald );
194  lcmt00 = lcmt;
195  lcmt -= pmb;
196  mblks = mblkd;
197  mblkd--;
198  Aoffi = ioffd;
199  ioffd += mbloc;
200  }
201 /*
202 * Scale the lower triangular part of sub( A ) when necessary.
203 */
204  tmp1 = m1 - ioffd + Aii - 1;
205  if( lower && ( tmp1 > 0 ) )
206  scal( C2F_CHAR( ALL ), &tmp1, &inbloc, &izero, ALPHA,
207  Mptr( A, ioffd+1, Aoffj+1, Ald, size ), &Ald );
208 
209  tmp1 = Aoffi - Aii + 1;
210  m1 -= tmp1;
211  n1 -= inbloc;
212  lcmt00 += low - ilow + qnb;
213  nblks--;
214  Aoffj += inbloc;
215 /*
216 * When the upper triangular part of sub( A ) should be scaled, take care of the
217 * n1 remaining columns of these tmp1 rows immediately.
218 */
219  if( upper && ( tmp1 > 0 ) && ( n1 > 0 ) )
220  scal( C2F_CHAR( ALL ), &tmp1, &n1, &izero, ALPHA,
221  Mptr( A, Aii, Aoffj+1, Ald, size ), &Ald );
222  Aii = Aoffi + 1;
223  Ajj = Aoffj + 1;
224  }
225  else if( GoEast )
226  {
227 /*
228 * Go one step east in the LCM table. Adjust the current LCM value as well as
229 * the local column index in A.
230 */
231  lcmt00 += low - ilow + qnb; nblks--; Aoffj += inbloc;
232 /*
233 * While there are blocks remaining that own lower entries, keep going east.
234 * Adjust the current LCM value as well as the local column index in A.
235 */
236  while( ( nblks > 0 ) && ( lcmt00 < low ) )
237  { lcmt00 += qnb; nblks--; Aoffj += Anb; }
238 /*
239 * Scale the lower triangular part of sub( A ) we just skipped when necessary.
240 */
241  tmp1 = MIN( Aoffj, jjmax ) - Ajj + 1;
242  if( lower && ( tmp1 > 0 ) )
243  {
244  scal( C2F_CHAR( ALL ), &m1, &tmp1, &izero, ALPHA,
245  Mptr( A, Aii, Ajj, Ald, size ), &Ald );
246  Ajj += tmp1;
247  n1 -= tmp1;
248  }
249 /*
250 * Return if no more column in the LCM table.
251 */
252  if( nblks <= 0 ) return;
253 /*
254 * lcmt00 >= low. The current block owns either diagonals or upper entries.
255 * Save the current position in the LCM table. After this row has been
256 * completely taken care of, re-start from this column and the next row of
257 * the LCM table.
258 */
259  lcmt = lcmt00; nblkd = nblks; joffd = Aoffj;
260 
261  nbloc = Anb;
262  while( ( nblkd > 0 ) && ( lcmt <= iupp ) )
263  {
264 /*
265 * A block owning diagonals lcmt00 >= low && lcmt00 <= iupp has been found.
266 */
267  if( nblkd == 1 ) nbloc = lnbloc;
268  scal( C2F_CHAR( UPLO ), &imbloc, &nbloc, &lcmt, ALPHA,
269  Mptr( A, Aii, joffd+1, Ald, size ), &Ald );
270  lcmt00 = lcmt;
271  lcmt += qnb;
272  nblks = nblkd;
273  nblkd--;
274  Aoffj = joffd;
275  joffd += nbloc;
276  }
277 /*
278 * Scale the upper triangular part of sub( A ) when necessary.
279 */
280  tmp1 = n1 - joffd + Ajj - 1;
281  if( upper && ( tmp1 > 0 ) )
282  scal( C2F_CHAR( ALL ), &imbloc, &tmp1, &izero, ALPHA,
283  Mptr( A, Aii, joffd+1, Ald, size ), &Ald );
284 
285  tmp1 = Aoffj - Ajj + 1;
286  m1 -= imbloc;
287  n1 -= tmp1;
288  lcmt00 -= ( iupp - upp + pmb );
289  mblks--;
290  Aoffi += imbloc;
291 /*
292 * When the lower triangular part of sub( A ) should be scaled, take care of the
293 * m1 remaining rows of these tmp1 columns immediately.
294 */
295  if( lower && ( m1 > 0 ) && ( tmp1 > 0 ) )
296  scal( C2F_CHAR( ALL ), &m1, &tmp1, &izero, ALPHA,
297  Mptr( A, Aoffi+1, Ajj, Ald, size ), &Ald );
298  Aii = Aoffi + 1;
299  Ajj = Aoffj + 1;
300  }
301 /*
302 * Loop over the remaining columns of the LCM table.
303 */
304  nbloc = Anb;
305  while( nblks > 0 )
306  {
307  if( nblks == 1 ) nbloc = lnbloc;
308 /*
309 * While there are blocks remaining that own upper entries, keep going south.
310 * Adjust the current LCM value as well as the local row index in A.
311 */
312  while( ( mblks > 0 ) && ( lcmt00 > upp ) )
313  { lcmt00 -= pmb; mblks--; Aoffi += Amb; }
314 /*
315 * Scale the upper triangular part of sub( A ) we just skipped when necessary.
316 */
317  tmp1 = MIN( Aoffi, iimax ) - Aii + 1;
318  if( upper && ( tmp1 > 0 ) )
319  {
320  scal( C2F_CHAR( ALL ), &tmp1, &n1, &izero, ALPHA,
321  Mptr( A, Aii, Aoffj+1, Ald, size ), &Ald );
322  Aii += tmp1;
323  m1 -= tmp1;
324  }
325 /*
326 * Return if no more row in the LCM table.
327 */
328  if( mblks <= 0 ) return;
329 /*
330 * lcmt00 <= upp. The current block owns either diagonals or lower entries.
331 * Save the current position in the LCM table. After this column has been
332 * completely taken care of, re-start from this row and the next column of
333 * the LCM table.
334 */
335  lcmt = lcmt00; mblkd = mblks; ioffd = Aoffi;
336 
337  mbloc = Amb;
338  while( ( mblkd > 0 ) && ( lcmt >= low ) )
339  {
340 /*
341 * A block owning diagonals lcmt00 >= low && lcmt00 <= upp has been found.
342 */
343  if( mblkd == 1 ) mbloc = lmbloc;
344  scal( C2F_CHAR( UPLO ), &mbloc, &nbloc, &lcmt, ALPHA,
345  Mptr( A, ioffd+1, Aoffj+1, Ald, size ), &Ald );
346  lcmt00 = lcmt;
347  lcmt -= pmb;
348  mblks = mblkd;
349  mblkd--;
350  Aoffi = ioffd;
351  ioffd += mbloc;
352  }
353 /*
354 * Scale the lower triangular part of sub( A ) when necessary.
355 */
356  tmp1 = m1 - ioffd + Aii - 1;
357  if( lower && ( tmp1 > 0 ) )
358  scal( C2F_CHAR( ALL ), &tmp1, &nbloc, &izero, ALPHA,
359  Mptr( A, ioffd+1, Aoffj+1, Ald, size ), &Ald );
360 
361  tmp1 = MIN( Aoffi, iimax ) - Aii + 1;
362  m1 -= tmp1;
363  n1 -= nbloc;
364  lcmt00 += qnb;
365  nblks--;
366  Aoffj += nbloc;
367 /*
368 * When the upper triangular part of sub( A ) should be scaled, take care of the
369 * n1 remaining columns of these tmp1 rows immediately.
370 */
371  if( upper && ( tmp1 > 0 ) && ( n1 > 0 ) )
372  scal( C2F_CHAR( ALL ), &tmp1, &n1, &izero, ALPHA,
373  Mptr( A, Aii, Aoffj+1, Ald, size ), &Ald );
374  Aii = Aoffi + 1;
375  Ajj = Aoffj + 1;
376  }
377 /*
378 * End of PB_Cplasca2
379 */
380 }
CCONJG
#define CCONJG
Definition: PBblas.h:21
TYPE
#define TYPE
Definition: clamov.c:7
MB_
#define MB_
Definition: PBtools.h:43
NB_
#define NB_
Definition: PBtools.h:44
PB_Cainfog2l
void PB_Cainfog2l()
LLD_
#define LLD_
Definition: PBtools.h:47
CNOCONJG
#define CNOCONJG
Definition: PBblas.h:19
PB_Cbinfo
void PB_Cbinfo()
CALL
#define CALL
Definition: PBblas.h:24
CLOWER
#define CLOWER
Definition: PBblas.h:25
ALL
#define ALL
Definition: PBblas.h:50
MIN
#define MIN(a_, b_)
Definition: PBtools.h:76
C2F_CHAR
#define C2F_CHAR(a)
Definition: pblas.h:121
TZSCAL_T
F_VOID_FCT(* TZSCAL_T)()
Definition: pblas.h:291
Cblacs_gridinfo
void Cblacs_gridinfo()
PBTYP_T
Definition: pblas.h:325
Mupcase
#define Mupcase(C)
Definition: PBtools.h:83
CUPPER
#define CUPPER
Definition: PBblas.h:26
Mptr
#define Mptr(a_, i_, j_, lda_, siz_)
Definition: PBtools.h:132
CTXT_
#define CTXT_
Definition: PBtools.h:38
PB_Cplasca2
void PB_Cplasca2(PBTYP_T *TYPE, char *UPLO, char *CONJUG, int M, int N, char *ALPHA, char *A, int IA, int JA, int *DESCA)
Definition: PB_Cplasca2.c:24