ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdlange.f
Go to the documentation of this file.
1  DOUBLE PRECISION FUNCTION pdlange( NORM, M, N, A, IA, JA, DESCA,
2  $ WORK )
3  IMPLICIT NONE
4 *
5 * -- ScaLAPACK auxiliary routine (version 1.7) --
6 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
7 * and University of California, Berkeley.
8 * May 1, 1997
9 *
10 * .. Scalar Arguments ..
11  CHARACTER norm
12  INTEGER ia, ja, m, n
13 * ..
14 * .. Array Arguments ..
15  INTEGER desca( * )
16  DOUBLE PRECISION a( * ), work( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * PDLANGE returns the value of the one norm, or the Frobenius norm,
23 * or the infinity norm, or the element of largest absolute value of a
24 * distributed matrix sub( A ) = A(IA:IA+M-1, JA:JA+N-1).
25 *
26 * PDLANGE returns the value
27 *
28 * ( max(abs(A(i,j))), NORM = 'M' or 'm' with IA <= i <= IA+M-1,
29 * ( and JA <= j <= JA+N-1,
30 * (
31 * ( norm1( sub( A ) ), NORM = '1', 'O' or 'o'
32 * (
33 * ( normI( sub( A ) ), NORM = 'I' or 'i'
34 * (
35 * ( normF( sub( A ) ), NORM = 'F', 'f', 'E' or 'e'
36 *
37 * where norm1 denotes the one norm of a matrix (maximum column sum),
38 * normI denotes the infinity norm of a matrix (maximum row sum) and
39 * normF denotes the Frobenius norm of a matrix (square root of sum of
40 * squares). Note that max(abs(A(i,j))) is not a matrix norm.
41 *
42 * Notes
43 * =====
44 *
45 * Each global data object is described by an associated description
46 * vector. This vector stores the information required to establish
47 * the mapping between an object element and its corresponding process
48 * and memory location.
49 *
50 * Let A be a generic term for any 2D block cyclicly distributed array.
51 * Such a global array has an associated description vector DESCA.
52 * In the following comments, the character _ should be read as
53 * "of the global array".
54 *
55 * NOTATION STORED IN EXPLANATION
56 * --------------- -------------- --------------------------------------
57 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
58 * DTYPE_A = 1.
59 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
60 * the BLACS process grid A is distribu-
61 * ted over. The context itself is glo-
62 * bal, but the handle (the integer
63 * value) may vary.
64 * M_A (global) DESCA( M_ ) The number of rows in the global
65 * array A.
66 * N_A (global) DESCA( N_ ) The number of columns in the global
67 * array A.
68 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
69 * the rows of the array.
70 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
71 * the columns of the array.
72 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
73 * row of the array A is distributed.
74 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
75 * first column of the array A is
76 * distributed.
77 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
78 * array. LLD_A >= MAX(1,LOCr(M_A)).
79 *
80 * Let K be the number of rows or columns of a distributed matrix,
81 * and assume that its process grid has dimension p x q.
82 * LOCr( K ) denotes the number of elements of K that a process
83 * would receive if K were distributed over the p processes of its
84 * process column.
85 * Similarly, LOCc( K ) denotes the number of elements of K that a
86 * process would receive if K were distributed over the q processes of
87 * its process row.
88 * The values of LOCr() and LOCc() may be determined via a call to the
89 * ScaLAPACK tool function, NUMROC:
90 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
91 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
92 * An upper bound for these quantities may be computed by:
93 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
94 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
95 *
96 * Arguments
97 * =========
98 *
99 * NORM (global input) CHARACTER
100 * Specifies the value to be returned in PDLANGE as described
101 * above.
102 *
103 * M (global input) INTEGER
104 * The number of rows to be operated on i.e the number of rows
105 * of the distributed submatrix sub( A ). When M = 0, PDLANGE
106 * is set to zero. M >= 0.
107 *
108 * N (global input) INTEGER
109 * The number of columns to be operated on i.e the number of
110 * columns of the distributed submatrix sub( A ). When N = 0,
111 * PDLANGE is set to zero. N >= 0.
112 *
113 * A (local input) DOUBLE PRECISION pointer into the local memory
114 * to an array of dimension (LLD_A, LOCc(JA+N-1)) containing the
115 * local pieces of the distributed matrix sub( A ).
116 *
117 * IA (global input) INTEGER
118 * The row index in the global array A indicating the first
119 * row of sub( A ).
120 *
121 * JA (global input) INTEGER
122 * The column index in the global array A indicating the
123 * first column of sub( A ).
124 *
125 * DESCA (global and local input) INTEGER array of dimension DLEN_.
126 * The array descriptor for the distributed matrix A.
127 *
128 * WORK (local workspace) DOUBLE PRECISION array dimension (LWORK)
129 * LWORK >= 0 if NORM = 'M' or 'm' (not referenced),
130 * Nq0 if NORM = '1', 'O' or 'o',
131 * Mp0 if NORM = 'I' or 'i',
132 * 0 if NORM = 'F', 'f', 'E' or 'e' (not referenced),
133 * where
134 *
135 * IROFFA = MOD( IA-1, MB_A ), ICOFFA = MOD( JA-1, NB_A ),
136 * IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ),
137 * IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ),
138 * Mp0 = NUMROC( M+IROFFA, MB_A, MYROW, IAROW, NPROW ),
139 * Nq0 = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ),
140 *
141 * INDXG2P and NUMROC are ScaLAPACK tool functions; MYROW,
142 * MYCOL, NPROW and NPCOL can be determined by calling the
143 * subroutine BLACS_GRIDINFO.
144 *
145 * =====================================================================
146 *
147 * .. Parameters ..
148  INTEGER block_cyclic_2d, csrc_, ctxt_, dlen_, dtype_,
149  $ lld_, mb_, m_, nb_, n_, rsrc_
150  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
151  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
152  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
153  DOUBLE PRECISION one, zero
154  parameter( one = 1.0d+0, zero = 0.0d+0 )
155 * ..
156 * .. Local Scalars ..
157  INTEGER i, iacol, iarow, ictxt, ii, icoff, ioffa,
158  $ iroff, j, jj, lda, mp, mycol, myrow, npcol,
159  $ nprow, nq
160  DOUBLE PRECISION sum, value
161 * ..
162 * .. Local Arrays ..
163  DOUBLE PRECISION ssq( 2 ), colssq( 2 )
164 * ..
165 * .. External Subroutines ..
166  EXTERNAL blacs_gridinfo, dcombssq, dgebr2d,
167  $ dgebs2d, dgamx2d, dgsum2d, dlassq,
168  $ infog2l, pdtreecomb
169 * ..
170 * .. External Functions ..
171  LOGICAL lsame
172  INTEGER idamax, numroc
173  EXTERNAL lsame, idamax, numroc
174 * ..
175 * .. Intrinsic Functions ..
176  INTRINSIC abs, max, min, mod, sqrt
177 * ..
178 * .. Executable Statements ..
179 *
180 * Get grid parameters.
181 *
182  ictxt = desca( ctxt_ )
183  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
184 *
185  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, ii, jj,
186  $ iarow, iacol )
187  iroff = mod( ia-1, desca( mb_ ) )
188  icoff = mod( ja-1, desca( nb_ ) )
189  mp = numroc( m+iroff, desca( mb_ ), myrow, iarow, nprow )
190  nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
191  IF( myrow.EQ.iarow )
192  $ mp = mp - iroff
193  IF( mycol.EQ.iacol )
194  $ nq = nq - icoff
195  lda = desca( lld_ )
196 *
197  IF( min( m, n ).EQ.0 ) THEN
198 *
199  VALUE = zero
200 *
201 ************************************************************************
202 * max norm
203 *
204  ELSE IF( lsame( norm, 'M' ) ) THEN
205 *
206 * Find max(abs(A(i,j))).
207 *
208  VALUE = zero
209  IF( nq.GT.0 .AND. mp.GT.0 ) THEN
210  ioffa = (jj-1)*lda
211  DO 20 j = jj, jj+nq-1
212  DO 10 i = ii, mp+ii-1
213  VALUE = max( VALUE, abs( a( ioffa+i ) ) )
214  10 CONTINUE
215  ioffa = ioffa + lda
216  20 CONTINUE
217  END IF
218  CALL dgamx2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, i, j, -1,
219  $ 0, 0 )
220 *
221 ************************************************************************
222 * one norm
223 *
224  ELSE IF( lsame( norm, 'O' ) .OR. norm.EQ.'1' ) THEN
225 *
226 * Find norm1( sub( A ) ).
227 *
228  IF( nq.GT.0 ) THEN
229  ioffa = ( jj - 1 ) * lda
230  DO 40 j = jj, jj+nq-1
231  sum = zero
232  IF( mp.GT.0 ) THEN
233  DO 30 i = ii, mp+ii-1
234  sum = sum + abs( a( ioffa+i ) )
235  30 CONTINUE
236  END IF
237  ioffa = ioffa + lda
238  work( j-jj+1 ) = sum
239  40 CONTINUE
240  END IF
241 *
242 * Find sum of global matrix columns and store on row 0 of
243 * process grid
244 *
245  CALL dgsum2d( ictxt, 'Columnwise', ' ', 1, nq, work, 1,
246  $ 0, mycol )
247 *
248 * Find maximum sum of columns for 1-norm
249 *
250  IF( myrow.EQ.0 ) THEN
251  IF( nq.GT.0 ) THEN
252  VALUE = work( idamax( nq, work, 1 ) )
253  ELSE
254  VALUE = zero
255  END IF
256  CALL dgamx2d( ictxt, 'Rowwise', ' ', 1, 1, VALUE, 1, i, j,
257  $ -1, 0, 0 )
258  END IF
259 *
260 ************************************************************************
261 * inf norm
262 *
263  ELSE IF( lsame( norm, 'I' ) ) THEN
264 *
265 * Find normI( sub( A ) ).
266 *
267  IF( mp.GT.0 ) THEN
268  ioffa = ii + ( jj - 1 ) * lda
269  DO 60 i = ii, ii+mp-1
270  sum = zero
271  IF( nq.GT.0 ) THEN
272  DO 50 j = ioffa, ioffa + nq*lda - 1, lda
273  sum = sum + abs( a( j ) )
274  50 CONTINUE
275  END IF
276  work( i-ii+1 ) = sum
277  ioffa = ioffa + 1
278  60 CONTINUE
279  END IF
280 *
281 * Find sum of global matrix rows and store on column 0 of
282 * process grid
283 *
284  CALL dgsum2d( ictxt, 'Rowwise', ' ', mp, 1, work, max( 1, mp ),
285  $ myrow, 0 )
286 *
287 * Find maximum sum of rows for supnorm
288 *
289  IF( mycol.EQ.0 ) THEN
290  IF( mp.GT.0 ) THEN
291  VALUE = work( idamax( mp, work, 1 ) )
292  ELSE
293  VALUE = zero
294  END IF
295  CALL dgamx2d( ictxt, 'Columnwise', ' ', 1, 1, VALUE, 1, i,
296  $ j, -1, 0, 0 )
297  END IF
298 *
299 ************************************************************************
300 * Frobenius norm
301 * SSQ(1) is scale
302 * SSQ(2) is sum-of-squares
303 *
304  ELSE IF( ( lsame( norm, 'F' ) ) .OR. ( lsame( norm, 'E' ) ) ) THEN
305 *
306 * Find normF( sub( A ) ).
307 *
308  ssq(1) = zero
309  ssq(2) = one
310  ioffa = ii + ( jj - 1 ) * lda
311  IF( nq.GT.0 ) THEN
312  DO 70 j = ioffa, ioffa + nq*lda - 1, lda
313  colssq(1) = zero
314  colssq(2) = one
315  CALL dlassq( mp, a( j ), 1, colssq(1), colssq(2) )
316  CALL dcombssq( ssq, colssq )
317  70 CONTINUE
318  END IF
319 *
320 * Perform the global scaled sum
321 *
322  CALL pdtreecomb( ictxt, 'All', 2, ssq, 0, 0, dcombssq )
323  VALUE = ssq( 1 ) * sqrt( ssq( 2 ) )
324 *
325  END IF
326 *
327  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
328  CALL dgebs2d( ictxt, 'All', ' ', 1, 1, VALUE, 1 )
329  ELSE
330  CALL dgebr2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, 0, 0 )
331  END IF
332 *
333  pdlange = VALUE
334 *
335  RETURN
336 *
337 * End of PDLANGE
338 *
339  END
max
#define max(A, B)
Definition: pcgemr.c:180
pdtreecomb
subroutine pdtreecomb(ICTXT, SCOPE, N, MINE, RDEST0, CDEST0, SUBPTR)
Definition: pdtreecomb.f:3
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
lsame
logical function lsame(CA, CB)
Definition: tools.f:1724
pdlange
double precision function pdlange(NORM, M, N, A, IA, JA, DESCA, WORK)
Definition: pdlange.f:3
numroc
integer function numroc(N, NB, IPROC, ISRCPROC, NPROCS)
Definition: numroc.f:2
dcombssq
subroutine dcombssq(V1, V2)
Definition: pdtreecomb.f:259
min
#define min(A, B)
Definition: pcgemr.c:181