ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pzgeqpf.f
Go to the documentation of this file.
1  SUBROUTINE pzgeqpf( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
2  $ LWORK, RWORK, LRWORK, INFO )
3 *
4 * -- ScaLAPACK routine (version 2.1) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * November 20, 2019
8 *
9 * .. Scalar Arguments ..
10  INTEGER IA, JA, INFO, LRWORK, LWORK, M, N
11 * ..
12 * .. Array Arguments ..
13  INTEGER DESCA( * ), IPIV( * )
14  DOUBLE PRECISION RWORK( * )
15  COMPLEX*16 A( * ), TAU( * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * PZGEQPF computes a QR factorization with column pivoting of a
22 * M-by-N distributed matrix sub( A ) = A(IA:IA+M-1,JA:JA+N-1):
23 *
24 * sub( A ) * P = Q * R.
25 *
26 * Notes
27 * =====
28 *
29 * Each global data object is described by an associated description
30 * vector. This vector stores the information required to establish
31 * the mapping between an object element and its corresponding process
32 * and memory location.
33 *
34 * Let A be a generic term for any 2D block cyclicly distributed array.
35 * Such a global array has an associated description vector DESCA.
36 * In the following comments, the character _ should be read as
37 * "of the global array".
38 *
39 * NOTATION STORED IN EXPLANATION
40 * --------------- -------------- --------------------------------------
41 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
42 * DTYPE_A = 1.
43 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
44 * the BLACS process grid A is distribu-
45 * ted over. The context itself is glo-
46 * bal, but the handle (the integer
47 * value) may vary.
48 * M_A (global) DESCA( M_ ) The number of rows in the global
49 * array A.
50 * N_A (global) DESCA( N_ ) The number of columns in the global
51 * array A.
52 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
53 * the rows of the array.
54 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
55 * the columns of the array.
56 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
57 * row of the array A is distributed.
58 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
59 * first column of the array A is
60 * distributed.
61 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
62 * array. LLD_A >= MAX(1,LOCr(M_A)).
63 *
64 * Let K be the number of rows or columns of a distributed matrix,
65 * and assume that its process grid has dimension p x q.
66 * LOCr( K ) denotes the number of elements of K that a process
67 * would receive if K were distributed over the p processes of its
68 * process column.
69 * Similarly, LOCc( K ) denotes the number of elements of K that a
70 * process would receive if K were distributed over the q processes of
71 * its process row.
72 * The values of LOCr() and LOCc() may be determined via a call to the
73 * ScaLAPACK tool function, NUMROC:
74 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
75 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
76 * An upper bound for these quantities may be computed by:
77 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
78 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
79 *
80 * Arguments
81 * =========
82 *
83 * M (global input) INTEGER
84 * The number of rows to be operated on, i.e. the number of rows
85 * of the distributed submatrix sub( A ). M >= 0.
86 *
87 * N (global input) INTEGER
88 * The number of columns to be operated on, i.e. the number of
89 * columns of the distributed submatrix sub( A ). N >= 0.
90 *
91 * A (local input/local output) COMPLEX*16 pointer into the
92 * local memory to an array of dimension (LLD_A, LOCc(JA+N-1)).
93 * On entry, the local pieces of the M-by-N distributed matrix
94 * sub( A ) which is to be factored. On exit, the elements on
95 * and above the diagonal of sub( A ) contain the min(M,N) by N
96 * upper trapezoidal matrix R (R is upper triangular if M >= N);
97 * the elements below the diagonal, with the array TAU, repre-
98 * sent the unitary matrix Q as a product of elementary
99 * reflectors (see Further Details).
100 *
101 * IA (global input) INTEGER
102 * The row index in the global array A indicating the first
103 * row of sub( A ).
104 *
105 * JA (global input) INTEGER
106 * The column index in the global array A indicating the
107 * first column of sub( A ).
108 *
109 * DESCA (global and local input) INTEGER array of dimension DLEN_.
110 * The array descriptor for the distributed matrix A.
111 *
112 * IPIV (local output) INTEGER array, dimension LOCc(JA+N-1).
113 * On exit, if IPIV(I) = K, the local i-th column of sub( A )*P
114 * was the global K-th column of sub( A ). IPIV is tied to the
115 * distributed matrix A.
116 *
117 * TAU (local output) COMPLEX*16, array, dimension
118 * LOCc(JA+MIN(M,N)-1). This array contains the scalar factors
119 * TAU of the elementary reflectors. TAU is tied to the
120 * distributed matrix A.
121 *
122 * WORK (local workspace/local output) COMPLEX*16 array,
123 * dimension (LWORK)
124 * On exit, WORK(1) returns the minimal and optimal LWORK.
125 *
126 * LWORK (local or global input) INTEGER
127 * The dimension of the array WORK.
128 * LWORK is local input and must be at least
129 * LWORK >= MAX(3,Mp0 + Nq0).
130 *
131 * If LWORK = -1, then LWORK is global input and a workspace
132 * query is assumed; the routine only calculates the minimum
133 * and optimal size for all work arrays. Each of these
134 * values is returned in the first entry of the corresponding
135 * work array, and no error message is issued by PXERBLA.
136 *
137 * RWORK (local workspace/local output) DOUBLE PRECISION array,
138 * dimension (LRWORK)
139 * On exit, RWORK(1) returns the minimal and optimal LRWORK.
140 *
141 * LRWORK (local or global input) INTEGER
142 * The dimension of the array RWORK.
143 * LRWORK is local input and must be at least
144 * LRWORK >= LOCc(JA+N-1)+Nq0.
145 *
146 * IROFF = MOD( IA-1, MB_A ), ICOFF = MOD( JA-1, NB_A ),
147 * IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ),
148 * IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ),
149 * Mp0 = NUMROC( M+IROFF, MB_A, MYROW, IAROW, NPROW ),
150 * Nq0 = NUMROC( N+ICOFF, NB_A, MYCOL, IACOL, NPCOL ),
151 * LOCc(JA+N-1) = NUMROC( JA+N-1, NB_A, MYCOL, CSRC_A, NPCOL )
152 *
153 * and NUMROC, INDXG2P are ScaLAPACK tool functions;
154 * MYROW, MYCOL, NPROW and NPCOL can be determined by calling
155 * the subroutine BLACS_GRIDINFO.
156 *
157 * If LRWORK = -1, then LRWORK is global input and a workspace
158 * query is assumed; the routine only calculates the minimum
159 * and optimal size for all work arrays. Each of these
160 * values is returned in the first entry of the corresponding
161 * work array, and no error message is issued by PXERBLA.
162 *
163 *
164 * INFO (global output) INTEGER
165 * = 0: successful exit
166 * < 0: If the i-th argument is an array and the j-entry had
167 * an illegal value, then INFO = -(i*100+j), if the i-th
168 * argument is a scalar and had an illegal value, then
169 * INFO = -i.
170 *
171 * Further Details
172 * ===============
173 *
174 * The matrix Q is represented as a product of elementary reflectors
175 *
176 * Q = H(1) H(2) . . . H(n)
177 *
178 * Each H(i) has the form
179 *
180 * H = I - tau * v * v'
181 *
182 * where tau is a complex scalar, and v is a complex vector with
183 * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
184 * A(ia+i-1:ia+m-1,ja+i-1).
185 *
186 * The matrix P is represented in jpvt as follows: If
187 * jpvt(j) = i
188 * then the jth column of P is the ith canonical unit vector.
189 *
190 * References
191 * ==========
192 *
193 * For modifications introduced in Scalapack 2.1
194 * LAWN 295
195 * New robust ScaLAPACK routine for computing the QR factorization with column pivoting
196 * Zvonimir Bujanovic, Zlatko Drmac
197 * http://www.netlib.org/lapack/lawnspdf/lawn295.pdf
198 *
199 * =====================================================================
200 *
201 * .. Parameters ..
202  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
203  $ lld_, mb_, m_, nb_, n_, rsrc_
204  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
205  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
206  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
207  DOUBLE PRECISION ONE, ZERO
208  parameter( one = 1.0d+0, zero = 0.0d+0 )
209 * ..
210 * .. Local Scalars ..
211  LOGICAL LQUERY
212  INTEGER I, IACOL, IAROW, ICOFF, ICTXT, ICURROW,
213  $ icurcol, ii, iia, ioffa, ipcol, iroff, itemp,
214  $ j, jb, jj, jja, jjpvt, jn, kb, k, kk, kstart,
215  $ kstep, lda, ll, lrwmin, lwmin, mn, mp, mycol,
216  $ myrow, npcol, nprow, nq, nq0, pvt
217  DOUBLE PRECISION TEMP, TEMP2, TOL3Z
218  COMPLEX*16 AJJ, ALPHA
219 * ..
220 * .. Local Arrays ..
221  INTEGER DESCN( DLEN_ ), IDUM1( 2 ), IDUM2( 2 )
222 * ..
223 * .. External Subroutines ..
224  EXTERNAL blacs_gridinfo, chk1mat, descset, igerv2d,
225  $ igesd2d, infog1l, infog2l, pchk1mat, pdamax,
226  $ pdznrm2, pxerbla, pzelset,
227  $ pzlarfc, pzlarfg, zcopy, zgebr2d,
228  $ zgebs2d, zgerv2d, zgesd2d, zlarfg,
229  $ zswap
230 * ..
231 * .. External Functions ..
232  INTEGER ICEIL, INDXG2P, NUMROC
233  EXTERNAL iceil, indxg2p, numroc
234  DOUBLE PRECISION DLAMCH
235  EXTERNAL dlamch
236 * ..
237 * .. Intrinsic Functions ..
238  INTRINSIC abs, dcmplx, dconjg, idint, max, min, mod, sqrt
239 * ..
240 * .. Executable Statements ..
241 *
242 * Get grid parameters
243 *
244  ictxt = desca( ctxt_ )
245  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
246 *
247 * Test the input parameters
248 *
249  info = 0
250  IF( nprow.EQ.-1 ) THEN
251  info = -(600+ctxt_)
252  ELSE
253  CALL chk1mat( m, 1, n, 2, ia, ja, desca, 6, info )
254  IF( info.EQ.0 ) THEN
255  iroff = mod( ia-1, desca( mb_ ) )
256  icoff = mod( ja-1, desca( nb_ ) )
257  iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
258  $ nprow )
259  iacol = indxg2p( ja, desca( nb_ ), mycol, desca( csrc_ ),
260  $ npcol )
261  mp = numroc( m+iroff, desca( mb_ ), myrow, iarow, nprow )
262  nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
263  nq0 = numroc( ja+n-1, desca( nb_ ), mycol, desca( csrc_ ),
264  $ npcol )
265  lwmin = max( 3, mp + nq )
266  lrwmin = nq0 + nq
267 *
268  work( 1 ) = dcmplx( dble( lwmin ) )
269  rwork( 1 ) = dble( lrwmin )
270  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 )
271  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
272  info = -10
273  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
274  info = -12
275  END IF
276  END IF
277  IF( lwork.EQ.-1 ) THEN
278  idum1( 1 ) = -1
279  ELSE
280  idum1( 1 ) = 1
281  END IF
282  idum2( 1 ) = 10
283  IF( lrwork.EQ.-1 ) THEN
284  idum1( 2 ) = -1
285  ELSE
286  idum1( 2 ) = 1
287  END IF
288  idum2( 2 ) = 12
289  CALL pchk1mat( m, 1, n, 2, ia, ja, desca, 6, 2, idum1, idum2,
290  $ info )
291  END IF
292 *
293  IF( info.NE.0 ) THEN
294  CALL pxerbla( ictxt, 'PZGEQPF', -info )
295  RETURN
296  ELSE IF( lquery ) THEN
297  RETURN
298  END IF
299 *
300 * Quick return if possible
301 *
302  IF( m.EQ.0 .OR. n.EQ.0 )
303  $ RETURN
304 *
305  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
306  $ iarow, iacol )
307  IF( myrow.EQ.iarow )
308  $ mp = mp - iroff
309  IF( mycol.EQ.iacol )
310  $ nq = nq - icoff
311  mn = min( m, n )
312  tol3z = sqrt( dlamch('Epsilon') )
313 *
314 * Initialize the array of pivots
315 *
316  lda = desca( lld_ )
317  jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
318  kstep = npcol * desca( nb_ )
319 *
320  IF( mycol.EQ.iacol ) THEN
321 *
322 * Handle first block separately
323 *
324  jb = jn - ja + 1
325  DO 10 ll = jja, jja+jb-1
326  ipiv( ll ) = ja + ll - jja
327  10 CONTINUE
328  kstart = jn + kstep - desca( nb_ )
329 *
330 * Loop over remaining block of columns
331 *
332  DO 30 kk = jja+jb, jja+nq-1, desca( nb_ )
333  kb = min( jja+nq-kk, desca( nb_ ) )
334  DO 20 ll = kk, kk+kb-1
335  ipiv( ll ) = kstart+ll-kk+1
336  20 CONTINUE
337  kstart = kstart + kstep
338  30 CONTINUE
339  ELSE
340  kstart = jn + ( mod( mycol-iacol+npcol, npcol )-1 )*
341  $ desca( nb_ )
342  DO 50 kk = jja, jja+nq-1, desca( nb_ )
343  kb = min( jja+nq-kk, desca( nb_ ) )
344  DO 40 ll = kk, kk+kb-1
345  ipiv( ll ) = kstart+ll-kk+1
346  40 CONTINUE
347  kstart = kstart + kstep
348  50 CONTINUE
349  END IF
350 *
351 * Initialize partial column norms, handle first block separately
352 *
353  CALL descset( descn, 1, desca( n_ ), 1, desca( nb_ ), myrow,
354  $ desca( csrc_ ), ictxt, 1 )
355 *
356  jj = jja
357  IF( mycol.EQ.iacol ) THEN
358  DO 60 kk = 0, jb-1
359  CALL pdznrm2( m, rwork( jj+kk ), a, ia, ja+kk, desca, 1 )
360  rwork( nq+jj+kk ) = rwork( jj+kk )
361  60 CONTINUE
362  jj = jj + jb
363  END IF
364  icurcol = mod( iacol+1, npcol )
365 *
366 * Loop over the remaining blocks of columns
367 *
368  DO 80 j = jn+1, ja+n-1, desca( nb_ )
369  jb = min( ja+n-j, desca( nb_ ) )
370 *
371  IF( mycol.EQ.icurcol ) THEN
372  DO 70 kk = 0, jb-1
373  CALL pdznrm2( m, rwork( jj+kk ), a, ia, j+kk, desca, 1 )
374  rwork( nq+jj+kk ) = rwork( jj+kk )
375  70 CONTINUE
376  jj = jj + jb
377  END IF
378  icurcol = mod( icurcol+1, npcol )
379  80 CONTINUE
380 *
381 * Compute factorization
382 *
383  DO 120 j = ja, ja+mn-1
384  i = ia + j - ja
385 *
386  CALL infog1l( j, desca( nb_ ), npcol, mycol, desca( csrc_ ),
387  $ jj, icurcol )
388  k = ja + n - j
389  IF( k.GT.1 ) THEN
390  CALL pdamax( k, temp, pvt, rwork, 1, j, descn,
391  $ descn( m_ ) )
392  ELSE
393  pvt = j
394  END IF
395  IF( j.NE.pvt ) THEN
396  CALL infog1l( pvt, desca( nb_ ), npcol, mycol,
397  $ desca( csrc_ ), jjpvt, ipcol )
398  IF( icurcol.EQ.ipcol ) THEN
399  IF( mycol.EQ.icurcol ) THEN
400  CALL zswap( mp, a( iia+(jj-1)*lda ), 1,
401  $ a( iia+(jjpvt-1)*lda ), 1 )
402  itemp = ipiv( jjpvt )
403  ipiv( jjpvt ) = ipiv( jj )
404  ipiv( jj ) = itemp
405  rwork( jjpvt ) = rwork( jj )
406  rwork( nq+jjpvt ) = rwork( nq+jj )
407  END IF
408  ELSE
409  IF( mycol.EQ.icurcol ) THEN
410 *
411  CALL zgesd2d( ictxt, mp, 1, a( iia+(jj-1)*lda ), lda,
412  $ myrow, ipcol )
413  work( 1 ) = dcmplx( dble( ipiv( jj ) ) )
414  work( 2 ) = dcmplx( rwork( jj ) )
415  work( 3 ) = dcmplx( rwork( jj + nq ) )
416  CALL zgesd2d( ictxt, 3, 1, work, 3, myrow, ipcol )
417 *
418  CALL zgerv2d( ictxt, mp, 1, a( iia+(jj-1)*lda ), lda,
419  $ myrow, ipcol )
420  CALL igerv2d( ictxt, 1, 1, ipiv( jj ), 1, myrow,
421  $ ipcol )
422 *
423  ELSE IF( mycol.EQ.ipcol ) THEN
424 *
425  CALL zgesd2d( ictxt, mp, 1, a( iia+(jjpvt-1)*lda ),
426  $ lda, myrow, icurcol )
427  CALL igesd2d( ictxt, 1, 1, ipiv( jjpvt ), 1, myrow,
428  $ icurcol )
429 *
430  CALL zgerv2d( ictxt, mp, 1, a( iia+(jjpvt-1)*lda ),
431  $ lda, myrow, icurcol )
432  CALL zgerv2d( ictxt, 3, 1, work, 3, myrow, icurcol )
433  ipiv( jjpvt ) = idint( dble( work( 1 ) ) )
434  rwork( jjpvt ) = dble( work( 2 ) )
435  rwork( jjpvt+nq ) = dble( work( 3 ) )
436 *
437  END IF
438 *
439  END IF
440 *
441  END IF
442 *
443 * Generate elementary reflector H(i)
444 *
445  CALL infog1l( i, desca( mb_ ), nprow, myrow, desca( rsrc_ ),
446  $ ii, icurrow )
447  IF( desca( m_ ).EQ.1 ) THEN
448  IF( myrow.EQ.icurrow ) THEN
449  IF( mycol.EQ.icurcol ) THEN
450  ioffa = ii+(jj-1)*desca( lld_ )
451  ajj = a( ioffa )
452  CALL zlarfg( 1, ajj, a( ioffa ), 1, tau( jj ) )
453  IF( n.GT.1 ) THEN
454  alpha = cmplx( one ) - dconjg( tau( jj ) )
455  CALL zgebs2d( ictxt, 'Rowwise', ' ', 1, 1, alpha,
456  $ 1 )
457  CALL zscal( nq-jj, alpha, a( ioffa+desca( lld_ ) ),
458  $ desca( lld_ ) )
459  END IF
460  CALL zgebs2d( ictxt, 'Columnwise', ' ', 1, 1,
461  $ tau( jj ), 1 )
462  a( ioffa ) = ajj
463  ELSE
464  IF( n.GT.1 ) THEN
465  CALL zgebr2d( ictxt, 'Rowwise', ' ', 1, 1, alpha,
466  $ 1, icurrow, icurcol )
467  CALL zscal( nq-jj+1, alpha, a( i ), desca( lld_ ) )
468  END IF
469  END IF
470  ELSE IF( mycol.EQ.icurcol ) THEN
471  CALL zgebr2d( ictxt, 'Columnwise', ' ', 1, 1, tau( jj ),
472  $ 1, icurrow, icurcol )
473  END IF
474 *
475  ELSE
476 *
477  CALL pzlarfg( m-j+ja, ajj, i, j, a, min( i+1, ia+m-1 ), j,
478  $ desca, 1, tau )
479  IF( j.LT.ja+n-1 ) THEN
480 *
481 * Apply H(i) to A(ia+j-ja:ia+m-1,j+1:ja+n-1) from the left
482 *
483  CALL pzelset( a, i, j, desca, dcmplx( one ) )
484  CALL pzlarfc( 'Left', m-j+ja, ja+n-1-j, a, i, j, desca,
485  $ 1, tau, a, i, j+1, desca, work )
486  END IF
487  CALL pzelset( a, i, j, desca, ajj )
488 *
489  END IF
490 *
491 * Update partial columns norms
492 *
493  IF( mycol.EQ.icurcol )
494  $ jj = jj + 1
495  IF( mod( j, desca( nb_ ) ).EQ.0 )
496  $ icurcol = mod( icurcol+1, npcol )
497  IF( (jja+nq-jj).GT.0 ) THEN
498  IF( myrow.EQ.icurrow ) THEN
499  CALL zgebs2d( ictxt, 'Columnwise', ' ', 1, jja+nq-jj,
500  $ a( ii+( min( jja+nq-1, jj )-1 )*lda ),
501  $ lda )
502  CALL zcopy( jja+nq-jj, a( ii+( min( jja+nq-1, jj )
503  $ -1)*lda ), lda, work( min( jja+nq-1, jj ) ),
504  $ 1 )
505  ELSE
506  CALL zgebr2d( ictxt, 'Columnwise', ' ', jja+nq-jj, 1,
507  $ work( min( jja+nq-1, jj ) ), max( 1, nq ),
508  $ icurrow, mycol )
509  END IF
510  END IF
511 *
512  jn = min( iceil( j+1, desca( nb_ ) ) * desca( nb_ ),
513  $ ja + n - 1 )
514  IF( mycol.EQ.icurcol ) THEN
515  DO 90 ll = jj, jj + jn - j - 1
516  IF( rwork( ll ).NE.zero ) THEN
517  temp = abs( work( ll ) ) / rwork( ll )
518  temp = max( zero, ( one+temp )*( one-temp ) )
519  temp2 = temp * ( rwork( ll ) / rwork( nq+ll ) )**2
520  IF( temp2.LE.tol3z ) THEN
521  IF( ia+m-1.GT.i ) THEN
522  CALL pdznrm2( ia+m-i-1, rwork( ll ), a,
523  $ i+1, j+ll-jj+1, desca, 1 )
524  rwork( nq+ll ) = rwork( ll )
525  ELSE
526  rwork( ll ) = zero
527  rwork( nq+ll ) = zero
528  END IF
529  ELSE
530  rwork( ll ) = rwork( ll ) * sqrt( temp )
531  END IF
532  END IF
533  90 CONTINUE
534  jj = jj + jn - j
535  END IF
536  icurcol = mod( icurcol+1, npcol )
537 *
538  DO 110 k = jn+1, ja+n-1, desca( nb_ )
539  kb = min( ja+n-k, desca( nb_ ) )
540 *
541  IF( mycol.EQ.icurcol ) THEN
542  DO 100 ll = jj, jj+kb-1
543  IF( rwork(ll).NE.zero ) THEN
544  temp = abs( work( ll ) ) / rwork( ll )
545  temp = max( zero, ( one+temp )*( one-temp ) )
546  temp2 = temp * ( rwork( ll ) / rwork( nq+ll ) )**2
547  IF( temp2.LE.tol3z ) THEN
548  IF( ia+m-1.GT.i ) THEN
549  CALL pdznrm2( ia+m-i-1, rwork( ll ), a,
550  $ i+1, k+ll-jj, desca, 1 )
551  rwork( nq+ll ) = rwork( ll )
552  ELSE
553  rwork( ll ) = zero
554  rwork( nq+ll ) = zero
555  END IF
556  ELSE
557  rwork( ll ) = rwork( ll ) * sqrt( temp )
558  END IF
559  END IF
560  100 CONTINUE
561  jj = jj + kb
562  END IF
563  icurcol = mod( icurcol+1, npcol )
564 *
565  110 CONTINUE
566 *
567  120 CONTINUE
568 *
569  work( 1 ) = dcmplx( dble( lwmin ) )
570  rwork( 1 ) = dble( lrwmin )
571 *
572  RETURN
573 *
574 * End of PZGEQPF
575 *
576  END
cmplx
float cmplx[2]
Definition: pblas.h:132
max
#define max(A, B)
Definition: pcgemr.c:180
infog1l
subroutine infog1l(GINDX, NB, NPROCS, MYROC, ISRCPROC, LINDX, ROCSRC)
Definition: infog1l.f:3
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pzlarfg
subroutine pzlarfg(N, ALPHA, IAX, JAX, X, IX, JX, DESCX, INCX, TAU)
Definition: pzlarfg.f:3
pzlarfc
subroutine pzlarfc(SIDE, M, N, V, IV, JV, DESCV, INCV, TAU, C, IC, JC, DESCC, WORK)
Definition: pzlarfc.f:3
pchk1mat
subroutine pchk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:3
descset
subroutine descset(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD)
Definition: descset.f:3
pzgeqpf
subroutine pzgeqpf(M, N, A, IA, JA, DESCA, IPIV, TAU, WORK, LWORK, RWORK, LRWORK, INFO)
Definition: pzgeqpf.f:3
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
pzelset
subroutine pzelset(A, IA, JA, DESCA, ALPHA)
Definition: pzelset.f:2
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
min
#define min(A, B)
Definition: pcgemr.c:181