Skip to content

Commit 33c15e1

Browse files
julielangoujulielangou
committed
New robust ScaLAPACK routine for computing the QR factorization with column pivoting
Contribution: Zvonimir Bujanovic, Zlatko Drmac LAWN 295 New robust ScaLAPACK routine for computing the QR factorization with column pivoting Zvonimir Bujanovic, Zlatko Drmac http://www.netlib.org/lapack/lawnspdf/lawn295.pdf
1 parent cf8e3d8 commit 33c15e1

File tree

4 files changed

+100
-61
lines changed

4 files changed

+100
-61
lines changed

‎SRC/pcgeqpf.f‎

Lines changed: 25 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
SUBROUTINE PCGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
2-
$ LWORK, RWORK, LRWORK, INFO )
2+
$ LWORK, RWORK, LRWORK, INFO )
33
*
4-
* -- ScaLAPACK routine (version 1.7) --
4+
* -- ScaLAPACK routine (version 2.1) --
55
* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
66
* and University of California, Berkeley.
7-
* March 14, 2000
7+
* November 20, 2019
88
*
99
* .. Scalar Arguments ..
1010
INTEGER IA, JA, INFO, LRWORK, LWORK, M, N
@@ -187,6 +187,15 @@ SUBROUTINE PCGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
187187
* jpvt(j) = i
188188
* then the jth column of P is the ith canonical unit vector.
189189
*
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+
*
190199
* =====================================================================
191200
*
192201
* .. Parameters ..
@@ -205,7 +214,7 @@ SUBROUTINE PCGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
205214
$ J, JB, JJ, JJA, JJPVT, JN, KB, K, KK, KSTART,
206215
$ KSTEP, LDA, LL, LRWMIN, LWMIN, MN, MP, MYCOL,
207216
$ MYROW, NPCOL, NPROW, NQ, NQ0, PVT
208-
REAL TEMP, TEMP2
217+
REAL TEMP, TEMP2, TOL3Z
209218
COMPLEX AJJ, ALPHA
210219
* ..
211220
* .. Local Arrays ..
@@ -221,6 +230,8 @@ SUBROUTINE PCGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
221230
* .. External Functions ..
222231
INTEGER ICEIL, INDXG2P, NUMROC
223232
EXTERNAL ICEIL, INDXG2P, NUMROC
233+
REAL SLAMCH
234+
EXTERNAL SLAMCH
224235
* ..
225236
* .. Intrinsic Functions ..
226237
INTRINSIC ABS, CMPLX, CONJG, IFIX, MAX, MIN, MOD, SQRT
@@ -297,6 +308,7 @@ SUBROUTINE PCGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
297308
IF( MYCOL.EQ.IACOL )
298309
$ NQ = NQ - ICOFF
299310
MN = MIN( M, N )
311+
TOL3Z = SQRT( SLAMCH('Epsilon') )
300312
*
301313
* Initialize the array of pivots
302314
*
@@ -501,14 +513,13 @@ SUBROUTINE PCGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
501513
IF( MYCOL.EQ.ICURCOL ) THEN
502514
DO 90 LL = JJ, JJ + JN - J - 1
503515
IF( RWORK( LL ).NE.ZERO ) THEN
504-
TEMP = ONE-( ABS( WORK( LL ) ) / RWORK( LL ) )**2
505-
TEMP = MAX( TEMP, ZERO )
506-
TEMP2 = ONE + 0.05E+0*TEMP*
507-
$ ( RWORK( LL ) / RWORK( NQ+LL ) )**2
508-
IF( TEMP2.EQ.ONE ) THEN
516+
TEMP = ABS( WORK( LL ) ) / RWORK( LL )
517+
TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
518+
TEMP2 = TEMP * ( RWORK( LL ) / RWORK( NQ+LL ) )**2
519+
IF( TEMP2.LE.TOL3Z ) THEN
509520
IF( IA+M-1.GT.I ) THEN
510521
CALL PSCNRM2( IA+M-I-1, RWORK( LL ), A,
511-
$ I+1, J+LL-JJ, DESCA, 1 )
522+
$ I+1, J+LL-JJ+1, DESCA, 1 )
512523
RWORK( NQ+LL ) = RWORK( LL )
513524
ELSE
514525
RWORK( LL ) = ZERO
@@ -529,11 +540,10 @@ SUBROUTINE PCGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
529540
IF( MYCOL.EQ.ICURCOL ) THEN
530541
DO 100 LL = JJ, JJ+KB-1
531542
IF( RWORK(LL).NE.ZERO ) THEN
532-
TEMP = ONE-( ABS( WORK( LL ) ) / RWORK( LL ) )**2
533-
TEMP = MAX( TEMP, ZERO )
534-
TEMP2 = ONE + 0.05E+0*TEMP*
535-
$ ( RWORK( LL ) / RWORK( NQ+LL ) )**2
536-
IF( TEMP2.EQ.ONE ) THEN
543+
TEMP = ABS( WORK( LL ) ) / RWORK( LL )
544+
TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
545+
TEMP2 = TEMP * ( RWORK( LL ) / RWORK( NQ+LL ) )**2
546+
IF( TEMP2.LE.TOL3Z ) THEN
537547
IF( IA+M-1.GT.I ) THEN
538548
CALL PSCNRM2( IA+M-I-1, RWORK( LL ), A,
539549
$ I+1, K+LL-JJ, DESCA, 1 )

‎SRC/pdgeqpf.f‎

Lines changed: 25 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
SUBROUTINE PDGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
2-
$ LWORK, INFO )
2+
$ LWORK, INFO )
33
*
4-
* -- ScaLAPACK routine (version 1.7) --
4+
* -- ScaLAPACK routine (version 2.1) --
55
* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
66
* and University of California, Berkeley.
7-
* March 14, 2000
7+
* November 20, 2019
88
*
99
* .. Scalar Arguments ..
1010
INTEGER IA, JA, INFO, LWORK, M, N
@@ -170,6 +170,15 @@ SUBROUTINE PDGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
170170
* jpvt(j) = i
171171
* then the jth column of P is the ith canonical unit vector.
172172
*
173+
* References
174+
* ==========
175+
*
176+
* For modifications introduced in Scalapack 2.1
177+
* LAWN 295
178+
* New robust ScaLAPACK routine for computing the QR factorization with column pivoting
179+
* Zvonimir Bujanovic, Zlatko Drmac
180+
* http://www.netlib.org/lapack/lawnspdf/lawn295.pdf
181+
*
173182
* =====================================================================
174183
*
175184
* .. Parameters ..
@@ -188,7 +197,7 @@ SUBROUTINE PDGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
188197
$ IROFF, ITEMP, J, JB, JJ, JJA, JJPVT, JN, KB,
189198
$ K, KK, KSTART, KSTEP, LDA, LL, LWMIN, MN, MP,
190199
$ MYCOL, MYROW, NPCOL, NPROW, NQ, NQ0, PVT
191-
DOUBLE PRECISION AJJ, ALPHA, TEMP, TEMP2
200+
DOUBLE PRECISION AJJ, ALPHA, TEMP, TEMP2, TOL3Z
192201
* ..
193202
* .. Local Arrays ..
194203
INTEGER DESCN( DLEN_ ), IDUM1( 1 ), IDUM2( 1 )
@@ -204,6 +213,8 @@ SUBROUTINE PDGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
204213
* .. External Functions ..
205214
INTEGER ICEIL, INDXG2P, NUMROC
206215
EXTERNAL ICEIL, INDXG2P, NUMROC
216+
DOUBLE PRECISION DLAMCH
217+
EXTERNAL DLAMCH
207218
* ..
208219
* .. Intrinsic Functions ..
209220
INTRINSIC ABS, DBLE, IDINT, MAX, MIN, MOD, SQRT
@@ -269,6 +280,7 @@ SUBROUTINE PDGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
269280
IF( MYCOL.EQ.IACOL )
270281
$ NQ = NQ - ICOFF
271282
MN = MIN( M, N )
283+
TOL3Z = SQRT( DLAMCH('Epsilon') )
272284
*
273285
* Initialize the array of pivots
274286
*
@@ -477,12 +489,10 @@ SUBROUTINE PDGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
477489
IF( MYCOL.EQ.ICURCOL ) THEN
478490
DO 90 LL = JJ-1, JJ + JN - J - 2
479491
IF( WORK( IPN+LL ).NE.ZERO ) THEN
480-
TEMP = ONE-( ABS( WORK( IPW+LL ) ) /
481-
$ WORK( IPN+LL ) )**2
482-
TEMP = MAX( TEMP, ZERO )
483-
TEMP2 = ONE + 0.05D+0*TEMP*
484-
$ ( WORK( IPN+LL ) / WORK( IPN+NQ+LL ) )**2
485-
IF( TEMP2.EQ.ONE ) THEN
492+
TEMP = ABS( WORK( IPW+LL ) ) / WORK( IPN+LL )
493+
TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
494+
TEMP2 = TEMP*( WORK( IPN+LL ) / WORK( IPN+NQ+LL ) )**2
495+
IF( TEMP2.LE.TOL3Z ) THEN
486496
IF( IA+M-1.GT.I ) THEN
487497
CALL PDNRM2( IA+M-I-1, WORK( IPN+LL ), A, I+1,
488498
$ J+LL-JJ+2, DESCA, 1 )
@@ -506,12 +516,11 @@ SUBROUTINE PDGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
506516
IF( MYCOL.EQ.ICURCOL ) THEN
507517
DO 100 LL = JJ-1, JJ+KB-2
508518
IF( WORK( IPN+LL ).NE.ZERO ) THEN
509-
TEMP = ONE-( ABS( WORK( IPW+LL ) ) /
510-
$ WORK( IPN+LL ) )**2
511-
TEMP = MAX( TEMP, ZERO )
512-
TEMP2 = ONE + 0.05D+0*TEMP*
513-
$ ( WORK( IPN+LL ) / WORK( IPN+NQ+LL ) )**2
514-
IF( TEMP2.EQ.ONE ) THEN
519+
TEMP = ABS( WORK( IPW+LL ) ) / WORK( IPN+LL )
520+
TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
521+
TEMP2 = TEMP*
522+
$ ( WORK( IPN+LL ) / WORK( IPN+NQ+LL ) )**2
523+
IF( TEMP2.LE.TOL3Z ) THEN
515524
IF( IA+M-1.GT.I ) THEN
516525
CALL PDNRM2( IA+M-I-1, WORK( IPN+LL ), A,
517526
$ I+1, K+LL-JJ+1, DESCA, 1 )

‎SRC/psgeqpf.f‎

Lines changed: 25 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
SUBROUTINE PSGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
22
$ LWORK, INFO )
33
*
4-
* -- ScaLAPACK routine (version 1.7) --
4+
* -- ScaLAPACK routine (version 2.1) --
55
* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
66
* and University of California, Berkeley.
7-
* March 14, 2000
7+
* November 20, 2019
88
*
99
* .. Scalar Arguments ..
1010
INTEGER IA, JA, INFO, LWORK, M, N
@@ -170,8 +170,18 @@ SUBROUTINE PSGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
170170
* jpvt(j) = i
171171
* then the jth column of P is the ith canonical unit vector.
172172
*
173+
* References
174+
* ==========
175+
*
176+
* For modifications introduced in Scalapack 2.1
177+
* LAWN 295
178+
* New robust ScaLAPACK routine for computing the QR factorization with column pivoting
179+
* Zvonimir Bujanovic, Zlatko Drmac
180+
* http://www.netlib.org/lapack/lawnspdf/lawn295.pdf
181+
*
173182
* =====================================================================
174183
*
184+
*
175185
* .. Parameters ..
176186
INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
177187
$ LLD_, MB_, M_, NB_, N_, RSRC_
@@ -188,7 +198,7 @@ SUBROUTINE PSGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
188198
$ IROFF, ITEMP, J, JB, JJ, JJA, JJPVT, JN, KB,
189199
$ K, KK, KSTART, KSTEP, LDA, LL, LWMIN, MN, MP,
190200
$ MYCOL, MYROW, NPCOL, NPROW, NQ, NQ0, PVT
191-
REAL AJJ, ALPHA, TEMP, TEMP2
201+
REAL AJJ, ALPHA, TEMP, TEMP2, TOL3Z
192202
* ..
193203
* .. Local Arrays ..
194204
INTEGER DESCN( DLEN_ ), IDUM1( 1 ), IDUM2( 1 )
@@ -203,6 +213,8 @@ SUBROUTINE PSGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
203213
* .. External Functions ..
204214
INTEGER ICEIL, INDXG2P, NUMROC
205215
EXTERNAL ICEIL, INDXG2P, NUMROC
216+
REAL SLAMCH
217+
EXTERNAL SLAMCH
206218
* ..
207219
* .. Intrinsic Functions ..
208220
INTRINSIC ABS, IFIX, MAX, MIN, MOD, REAL, SQRT
@@ -268,6 +280,7 @@ SUBROUTINE PSGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
268280
IF( MYCOL.EQ.IACOL )
269281
$ NQ = NQ - ICOFF
270282
MN = MIN( M, N )
283+
TOL3Z = SQRT( SLAMCH('Epsilon') )
271284
*
272285
* Initialize the array of pivots
273286
*
@@ -476,12 +489,10 @@ SUBROUTINE PSGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
476489
IF( MYCOL.EQ.ICURCOL ) THEN
477490
DO 90 LL = JJ-1, JJ + JN - J - 2
478491
IF( WORK( IPN+LL ).NE.ZERO ) THEN
479-
TEMP = ONE-( ABS( WORK( IPW+LL ) ) /
480-
$ WORK( IPN+LL ) )**2
481-
TEMP = MAX( TEMP, ZERO )
482-
TEMP2 = ONE + 0.05E+0*TEMP*
483-
$ ( WORK( IPN+LL ) / WORK( IPN+NQ+LL ) )**2
484-
IF( TEMP2.EQ.ONE ) THEN
492+
TEMP = ABS( WORK( IPW+LL ) ) / WORK( IPN+LL )
493+
TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
494+
TEMP2 = TEMP*( WORK( IPN+LL ) / WORK( IPN+NQ+LL ) )**2
495+
IF( TEMP2.LE.TOL3Z ) THEN
485496
IF( IA+M-1.GT.I ) THEN
486497
CALL PSNRM2( IA+M-I-1, WORK( IPN+LL ), A, I+1,
487498
$ J+LL-JJ+2, DESCA, 1 )
@@ -505,12 +516,11 @@ SUBROUTINE PSGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
505516
IF( MYCOL.EQ.ICURCOL ) THEN
506517
DO 100 LL = JJ-1, JJ+KB-2
507518
IF( WORK( IPN+LL ).NE.ZERO ) THEN
508-
TEMP = ONE-( ABS( WORK( IPW+LL ) ) /
509-
$ WORK( IPN+LL ) )**2
510-
TEMP = MAX( TEMP, ZERO )
511-
TEMP2 = ONE + 0.05E+0*TEMP*
512-
$ ( WORK( IPN+LL ) / WORK( IPN+NQ+LL ) )**2
513-
IF( TEMP2.EQ.ONE ) THEN
519+
TEMP = ABS( WORK( IPW+LL ) ) / WORK( IPN+LL )
520+
TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
521+
TEMP2 = TEMP*
522+
$ ( WORK( IPN+LL ) / WORK( IPN+NQ+LL ) )**2
523+
IF( TEMP2.LE.TOL3Z ) THEN
514524
IF( IA+M-1.GT.I ) THEN
515525
CALL PSNRM2( IA+M-I-1, WORK( IPN+LL ), A,
516526
$ I+1, K+LL-JJ+1, DESCA, 1 )

‎SRC/pzgeqpf.f‎

Lines changed: 25 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
SUBROUTINE PZGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
2-
$ LWORK, RWORK, LRWORK, INFO )
2+
$ LWORK, RWORK, LRWORK, INFO )
33
*
4-
* -- ScaLAPACK routine (version 1.7) --
4+
* -- ScaLAPACK routine (version 2.1) --
55
* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
66
* and University of California, Berkeley.
7-
* March 14, 2000
7+
* November 20, 2019
88
*
99
* .. Scalar Arguments ..
1010
INTEGER IA, JA, INFO, LRWORK, LWORK, M, N
@@ -187,6 +187,15 @@ SUBROUTINE PZGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
187187
* jpvt(j) = i
188188
* then the jth column of P is the ith canonical unit vector.
189189
*
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+
*
190199
* =====================================================================
191200
*
192201
* .. Parameters ..
@@ -205,7 +214,7 @@ SUBROUTINE PZGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
205214
$ J, JB, JJ, JJA, JJPVT, JN, KB, K, KK, KSTART,
206215
$ KSTEP, LDA, LL, LRWMIN, LWMIN, MN, MP, MYCOL,
207216
$ MYROW, NPCOL, NPROW, NQ, NQ0, PVT
208-
DOUBLE PRECISION TEMP, TEMP2
217+
DOUBLE PRECISION TEMP, TEMP2, TOL3Z
209218
COMPLEX*16 AJJ, ALPHA
210219
* ..
211220
* .. Local Arrays ..
@@ -222,6 +231,8 @@ SUBROUTINE PZGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
222231
* .. External Functions ..
223232
INTEGER ICEIL, INDXG2P, NUMROC
224233
EXTERNAL ICEIL, INDXG2P, NUMROC
234+
DOUBLE PRECISION DLAMCH
235+
EXTERNAL DLAMCH
225236
* ..
226237
* .. Intrinsic Functions ..
227238
INTRINSIC ABS, DCMPLX, DCONJG, IDINT, MAX, MIN, MOD, SQRT
@@ -298,6 +309,7 @@ SUBROUTINE PZGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
298309
IF( MYCOL.EQ.IACOL )
299310
$ NQ = NQ - ICOFF
300311
MN = MIN( M, N )
312+
TOL3Z = SQRT( DLAMCH('Epsilon') )
301313
*
302314
* Initialize the array of pivots
303315
*
@@ -502,14 +514,13 @@ SUBROUTINE PZGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
502514
IF( MYCOL.EQ.ICURCOL ) THEN
503515
DO 90 LL = JJ, JJ + JN - J - 1
504516
IF( RWORK( LL ).NE.ZERO ) THEN
505-
TEMP = ONE-( ABS( WORK( LL ) ) / RWORK( LL ) )**2
506-
TEMP = MAX( TEMP, ZERO )
507-
TEMP2 = ONE + 0.05D+0*TEMP*
508-
$ ( RWORK( LL ) / RWORK( NQ+LL ) )**2
509-
IF( TEMP2.EQ.ONE ) 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
510521
IF( IA+M-1.GT.I ) THEN
511522
CALL PDZNRM2( IA+M-I-1, RWORK( LL ), A,
512-
$ I+1, J+LL-JJ, DESCA, 1 )
523+
$ I+1, J+LL-JJ+1, DESCA, 1 )
513524
RWORK( NQ+LL ) = RWORK( LL )
514525
ELSE
515526
RWORK( LL ) = ZERO
@@ -530,11 +541,10 @@ SUBROUTINE PZGEQPF( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
530541
IF( MYCOL.EQ.ICURCOL ) THEN
531542
DO 100 LL = JJ, JJ+KB-1
532543
IF( RWORK(LL).NE.ZERO ) THEN
533-
TEMP = ONE-( ABS( WORK( LL ) ) / RWORK( LL ) )**2
534-
TEMP = MAX( TEMP, ZERO )
535-
TEMP2 = ONE + 0.05D+0*TEMP*
536-
$ ( RWORK( LL ) / RWORK( NQ+LL ) )**2
537-
IF( TEMP2.EQ.ONE ) 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
538548
IF( IA+M-1.GT.I ) THEN
539549
CALL PDZNRM2( IA+M-I-1, RWORK( LL ), A,
540550
$ I+1, K+LL-JJ, DESCA, 1 )

0 commit comments

Comments
 (0)