1 /* ---------------------------------------------------------------------
2 *
3 * -- PBLAS 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__
pzaxpy_(int * N,double * ALPHA,double * X,int * IX,int * JX,int * DESCX,int * INCX,double * Y,int * IY,int * JY,int * DESCY,int * INCY)20 void pzaxpy_( int * N,
21 double * ALPHA,
22 double * X, int * IX, int * JX, int * DESCX, int * INCX,
23 double * Y, int * IY, int * JY, int * DESCY, int * INCY )
24 #else
25 void pzaxpy_( N, ALPHA, X, IX, JX, DESCX, INCX, Y, IY, JY, DESCY, INCY )
26 /*
27 * .. Scalar Arguments ..
28 */
29 int * INCX, * INCY, * IX, * IY, * JX, * JY, * N;
30 double * ALPHA;
31 /*
32 * .. Array Arguments ..
33 */
34 int * DESCX, * DESCY;
35 double * X, * Y;
36 #endif
37 {
38 /*
39 * Purpose
40 * =======
41 *
42 * PZAXPY adds one subvector to another,
43 *
44 * sub( Y ) := sub( Y ) + alpha * sub( X ),
45 *
46 * where
47 *
48 * sub( X ) denotes X(IX,JX:JX+N-1) if INCX = M_X,
49 * X(IX:IX+N-1,JX) if INCX = 1 and INCX <> M_X, and,
50 *
51 * sub( Y ) denotes Y(IY,JY:JY+N-1) if INCY = M_Y,
52 * Y(IY:IY+N-1,JY) if INCY = 1 and INCY <> M_Y.
53 *
54 * Notes
55 * =====
56 *
57 * A description vector is associated with each 2D block-cyclicly dis-
58 * tributed matrix. This vector stores the information required to
59 * establish the mapping between a matrix entry and its corresponding
60 * process and memory location.
61 *
62 * In the following comments, the character _ should be read as
63 * "of the distributed matrix". Let A be a generic term for any 2D
64 * block cyclicly distributed matrix. Its description vector is DESC_A:
65 *
66 * NOTATION STORED IN EXPLANATION
67 * ---------------- --------------- ------------------------------------
68 * DTYPE_A (global) DESCA[ DTYPE_ ] The descriptor type.
69 * CTXT_A (global) DESCA[ CTXT_ ] The BLACS context handle, indicating
70 * the NPROW x NPCOL BLACS process grid
71 * A is distributed over. The context
72 * itself is global, but the handle
73 * (the integer value) may vary.
74 * M_A (global) DESCA[ M_ ] The number of rows in the distribu-
75 * ted matrix A, M_A >= 0.
76 * N_A (global) DESCA[ N_ ] The number of columns in the distri-
77 * buted matrix A, N_A >= 0.
78 * IMB_A (global) DESCA[ IMB_ ] The number of rows of the upper left
79 * block of the matrix A, IMB_A > 0.
80 * INB_A (global) DESCA[ INB_ ] The number of columns of the upper
81 * left block of the matrix A,
82 * INB_A > 0.
83 * MB_A (global) DESCA[ MB_ ] The blocking factor used to distri-
84 * bute the last M_A-IMB_A rows of A,
85 * MB_A > 0.
86 * NB_A (global) DESCA[ NB_ ] The blocking factor used to distri-
87 * bute the last N_A-INB_A columns of
88 * A, NB_A > 0.
89 * RSRC_A (global) DESCA[ RSRC_ ] The process row over which the first
90 * row of the matrix A is distributed,
91 * NPROW > RSRC_A >= 0.
92 * CSRC_A (global) DESCA[ CSRC_ ] The process column over which the
93 * first column of A is distributed.
94 * NPCOL > CSRC_A >= 0.
95 * LLD_A (local) DESCA[ LLD_ ] The leading dimension of the local
96 * array storing the local blocks of
97 * the distributed matrix A,
98 * IF( Lc( 1, N_A ) > 0 )
99 * LLD_A >= MAX( 1, Lr( 1, M_A ) )
100 * ELSE
101 * LLD_A >= 1.
102 *
103 * Let K be the number of rows of a matrix A starting at the global in-
104 * dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
105 * that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
106 * receive if these K rows were distributed over NPROW processes. If K
107 * is the number of columns of a matrix A starting at the global index
108 * JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number of co-
109 * lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would receive if
110 * these K columns were distributed over NPCOL processes.
111 *
112 * The values of Lr() and Lc() may be determined via a call to the func-
113 * tion PB_Cnumroc:
114 * Lr( IA, K ) = PB_Cnumroc( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
115 * Lc( JA, K ) = PB_Cnumroc( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
116 *
117 * Arguments
118 * =========
119 *
120 * N (global input) INTEGER.
121 * On entry, N specifies the length of the subvectors to be
122 * added. N must be at least zero.
123 *
124 * ALPHA (global input) COMPLEX*16
125 * On entry, ALPHA specifies the scalar alpha. When ALPHA is
126 * supplied as zero then the local entries of the array X cor-
127 * responding to the entries of the subvector sub( X ) need not
128 * be set on input.
129 *
130 * X (local input) COMPLEX*16 array
131 * On entry, X is an array of dimension (LLD_X, Kx), where LLD_X
132 * is at least MAX( 1, Lr( 1, IX ) ) when INCX = M_X and
133 * MAX( 1, Lr( 1, IX+N-1 ) ) otherwise, and, Kx is at least
134 * Lc( 1, JX+N-1 ) when INCX = M_X and Lc( 1, JX ) otherwise.
135 * Before entry, this array contains the local entries of the
136 * matrix X.
137 *
138 * IX (global input) INTEGER
139 * On entry, IX specifies X's global row index, which points to
140 * the beginning of the submatrix sub( X ).
141 *
142 * JX (global input) INTEGER
143 * On entry, JX specifies X's global column index, which points
144 * to the beginning of the submatrix sub( X ).
145 *
146 * DESCX (global and local input) INTEGER array
147 * On entry, DESCX is an integer array of dimension DLEN_. This
148 * is the array descriptor for the matrix X.
149 *
150 * INCX (global input) INTEGER
151 * On entry, INCX specifies the global increment for the
152 * elements of X. Only two values of INCX are supported in
153 * this version, namely 1 and M_X. INCX must not be zero.
154 *
155 * Y (local input/local output) COMPLEX*16 array
156 * On entry, Y is an array of dimension (LLD_Y, Ky), where LLD_Y
157 * is at least MAX( 1, Lr( 1, IY ) ) when INCY = M_Y and
158 * MAX( 1, Lr( 1, IY+N-1 ) ) otherwise, and, Ky is at least
159 * Lc( 1, JY+N-1 ) when INCY = M_Y and Lc( 1, JY ) otherwise.
160 * Before entry, this array contains the local entries of the
161 * matrix Y. On exit, sub( Y ) is overwritten with the updated
162 * subvector.
163 *
164 * IY (global input) INTEGER
165 * On entry, IY specifies Y's global row index, which points to
166 * the beginning of the submatrix sub( Y ).
167 *
168 * JY (global input) INTEGER
169 * On entry, JY specifies Y's global column index, which points
170 * to the beginning of the submatrix sub( Y ).
171 *
172 * DESCY (global and local input) INTEGER array
173 * On entry, DESCY is an integer array of dimension DLEN_. This
174 * is the array descriptor for the matrix Y.
175 *
176 * INCY (global input) INTEGER
177 * On entry, INCY specifies the global increment for the
178 * elements of Y. Only two values of INCY are supported in
179 * this version, namely 1 and M_Y. INCY must not be zero.
180 *
181 * -- Written on April 1, 1998 by
182 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
183 *
184 * ---------------------------------------------------------------------
185 */
186 /*
187 * .. Local Scalars ..
188 */
189 int Xi, Xj, Yi, Yj, ctxt, info, mycol, myrow, npcol, nprow;
190 PBTYP_T * type;
191 /*
192 * .. Local Arrays ..
193 */
194 int Xd[DLEN_], Yd[DLEN_];
195 /* ..
196 * .. Executable Statements ..
197 *
198 */
199 PB_CargFtoC( *IX, *JX, DESCX, &Xi, &Xj, Xd );
200 PB_CargFtoC( *IY, *JY, DESCY, &Yi, &Yj, Yd );
201 #ifndef NO_ARGCHK
202 /*
203 * Test the input parameters
204 */
205 Cblacs_gridinfo( ( ctxt = Xd[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
206 if( !( info = ( ( nprow == -1 ) ? -( 601 + CTXT_ ) : 0 ) ) )
207 {
208 PB_Cchkvec( ctxt, "PZAXPY", "X", *N, 1, Xi, Xj, Xd, *INCX, 6, &info );
209 PB_Cchkvec( ctxt, "PZAXPY", "Y", *N, 1, Yi, Yj, Yd, *INCY, 11, &info );
210 }
211 if( info ) { PB_Cabort( ctxt, "PZAXPY", info ); return; }
212 #endif
213 /*
214 * Quick return if possible
215 */
216 if( ( *N == 0 ) ||
217 ( ( ALPHA[REAL_PART] == ZERO ) && ( ALPHA[IMAG_PART] == ZERO ) ) )
218 return;
219 /*
220 * Get type structure
221 */
222 type = PB_Cztypeset();
223 /*
224 * Start the operations
225 */
226 if( *INCX == Xd[M_] )
227 {
228 PB_Cpaxpby( type, NOCONJG, 1, *N, ((char *) ALPHA), ((char *) X),
229 Xi, Xj, Xd, ROW, type->one, ((char *) Y), Yi, Yj, Yd,
230 ( *INCY == Yd[M_] ? ROW : COLUMN ) );
231 }
232 else
233 {
234 PB_Cpaxpby( type, NOCONJG, *N, 1, ((char *) ALPHA), ((char *) X),
235 Xi, Xj, Xd, COLUMN, type->one, ((char *) Y), Yi, Yj, Yd,
236 ( *INCY == Yd[M_] ? ROW : COLUMN ) );
237 }
238 /*
239 * End of PZAXPY
240 */
241 }
242