--- rpl/lapack/lapack/zgejsv.f 2015/11/26 11:44:21 1.1 +++ rpl/lapack/lapack/zgejsv.f 2016/08/27 15:27:12 1.2 @@ -27,7 +27,7 @@ * INTEGER INFO, LDA, LDU, LDV, LWORK, M, N * .. * .. Array Arguments .. -* DOUBLE COMPLEX A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( LWORK ) +* COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( LWORK ) * DOUBLE PRECISION SVA( N ), RWORK( LRWORK ) * INTEGER IWORK( * ) * CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV @@ -39,21 +39,22 @@ *> *> \verbatim *> -* ZGEJSV computes the singular value decomposition (SVD) of a real M-by-N -* matrix [A], where M >= N. The SVD of [A] is written as -* -* [A] = [U] * [SIGMA] * [V]^*, -* -* where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N -* diagonal elements, [U] is an M-by-N (or M-by-M) orthonormal matrix, and -* [V] is an N-by-N orthogonal matrix. The diagonal elements of [SIGMA] are -* the singular values of [A]. The columns of [U] and [V] are the left and -* the right singular vectors of [A], respectively. The matrices [U] and [V] -* are computed and stored in the arrays U and V, respectively. The diagonal -* of [SIGMA] is computed and stored in the array SVA. -* -* Arguments: -* ========== +*> ZGEJSV computes the singular value decomposition (SVD) of a complex M-by-N +*> matrix [A], where M >= N. The SVD of [A] is written as +*> +*> [A] = [U] * [SIGMA] * [V]^*, +*> +*> where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N +*> diagonal elements, [U] is an M-by-N (or M-by-M) unitary matrix, and +*> [V] is an N-by-N unitary matrix. The diagonal elements of [SIGMA] are +*> the singular values of [A]. The columns of [U] and [V] are the left and +*> the right singular vectors of [A], respectively. The matrices [U] and [V] +*> are computed and stored in the arrays U and V, respectively. The diagonal +*> of [SIGMA] is computed and stored in the array SVA. +*> \endverbatim +*> +*> Arguments: +*> ========== *> *> \param[in] JOBA *> \verbatim @@ -193,7 +194,7 @@ *> *> \param[in,out] A *> \verbatim -*> A is DOUBLE COMPLEX array, dimension (LDA,N) +*> A is COMPLEX*16 array, dimension (LDA,N) *> On entry, the M-by-N matrix A. *> \endverbatim *> @@ -221,7 +222,7 @@ *> *> \param[out] U *> \verbatim -*> U is DOUBLE COMPLEX array, dimension ( LDU, N ) +*> U is COMPLEX*16 array, dimension ( LDU, N ) *> If JOBU = 'U', then U contains on exit the M-by-N matrix of *> the left singular vectors. *> If JOBU = 'F', then U contains on exit the M-by-M matrix of @@ -234,7 +235,7 @@ *> copied back to the V array. This 'W' option is just *> a reminder to the caller that in this case U is *> reserved as workspace of length N*N. -*> If JOBU = 'N' U is not referenced. +*> If JOBU = 'N' U is not referenced, unless JOBT='T'. *> \endverbatim *> *> \param[in] LDU @@ -246,7 +247,7 @@ *> *> \param[out] V *> \verbatim -*> V is DOUBLE COMPLEX array, dimension ( LDV, N ) +*> V is COMPLEX*16 array, dimension ( LDV, N ) *> If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of *> the right singular vectors; *> If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N), @@ -256,7 +257,7 @@ *> copied back to the U array. This 'W' option is just *> a reminder to the caller that in this case V is *> reserved as workspace of length N*N. -*> If JOBV = 'N' V is not referenced. +*> If JOBV = 'N' V is not referenced, unless JOBT='T'. *> \endverbatim *> *> \param[in] LDV @@ -268,8 +269,7 @@ *> *> \param[out] CWORK *> \verbatim -*> CWORK (workspace) -*> CWORK is DOUBLE COMPLEX array, dimension at least LWORK. +*> CWORK is COMPLEX*16 array, dimension at least LWORK. *> \endverbatim *> *> \param[in] LWORK @@ -279,7 +279,7 @@ *> LWORK depends on the job: *> *> 1. If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and -*> 1.1 .. no scaled condition estimate required (JOBE.EQ.'N'): +*> 1.1 .. no scaled condition estimate required (JOBA.NE.'E'.AND.JOBA.NE.'G'): *> LWORK >= 2*N+1. This is the minimal requirement. *> ->> For optimal performance (blocked code) the optimal value *> is LWORK >= N + (N+1)*NB. Here NB is the optimal @@ -293,33 +293,33 @@ *> is LWORK >= max(N+(N+1)*NB, N*N+3*N). *> In general, the optimal length LWORK is computed as *> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), -*> N+N*N+LWORK(CPOCON)). +*> N+N*N+LWORK(ZPOCON)). *> *> 2. If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'), *> (JOBU.EQ.'N') *> -> the minimal requirement is LWORK >= 3*N. *> -> For optimal performance, LWORK >= max(N+(N+1)*NB, 3*N,2*N+N*NB), -*> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQ, -*> CUNMLQ. In general, the optimal length LWORK is computed as -*> LWORK >= max(N+LWORK(ZGEQP3), N+LWORK(CPOCON), N+LWORK(ZGESVJ), -*> N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(CUNMLQ)). +*> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQF, +*> ZUNMLQ. In general, the optimal length LWORK is computed as +*> LWORK >= max(N+LWORK(ZGEQP3), N+LWORK(ZPOCON), N+LWORK(ZGESVJ), +*> N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)). *> *> 3. If SIGMA and the left singular vectors are needed *> -> the minimal requirement is LWORK >= 3*N. *> -> For optimal performance: *> if JOBU.EQ.'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB), -*> where NB is the optimal block size for ZGEQP3, ZGEQRF, CUNMQR. +*> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR. *> In general, the optimal length LWORK is computed as -*> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(CPOCON), -*> 2*N+LWORK(ZGEQRF), N+LWORK(CUNMQR)). +*> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZPOCON), +*> 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)). *> *> 4. If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and *> 4.1. if JOBV.EQ.'V' *> the minimal requirement is LWORK >= 5*N+2*N*N. *> 4.2. if JOBV.EQ.'J' the minimal requirement is *> LWORK >= 4*N+N*N. -*> In both cases, the allocated CWORK can accomodate blocked runs -*> of ZGEQP3, ZGEQRF, ZGELQF, SUNMQR, CUNMLQ. +*> In both cases, the allocated CWORK can accommodate blocked runs +*> of ZGEQP3, ZGEQRF, ZGELQF, ZUNMQR, ZUNMLQ. *> \endverbatim *> *> \param[out] RWORK @@ -433,7 +433,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date November 2015 +*> \date June 2016 * *> \ingroup complex16GEsing * @@ -491,19 +491,19 @@ *> *> \verbatim *> -* [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. -* SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. -* LAPACK Working note 169. -* [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. -* SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. -* LAPACK Working note 170. -* [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR -* factorization software - a case study. -* ACM Trans. math. Softw. Vol. 35, No 2 (2008), pp. 1-28. -* LAPACK Working note 176. -* [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, -* QSVD, (H,K)-SVD computations. -* Department of Mathematics, University of Zagreb, 2008. +*> [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. +*> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. +*> LAPACK Working note 169. +*> [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. +*> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. +*> LAPACK Working note 170. +*> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR +*> factorization software - a case study. +*> ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28. +*> LAPACK Working note 176. +*> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, +*> QSVD, (H,K)-SVD computations. +*> Department of Mathematics, University of Zagreb, 2008. *> \endverbatim * *> \par Bugs, examples and comments: @@ -517,17 +517,17 @@ $ M, N, A, LDA, SVA, U, LDU, V, LDV, $ CWORK, LWORK, RWORK, LRWORK, IWORK, INFO ) * -* -- LAPACK computational routine (version 3.6.0) -- +* -- LAPACK computational routine (version 3.6.1) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2015 +* June 2016 * * .. Scalar Arguments .. IMPLICIT NONE INTEGER INFO, LDA, LDU, LDV, LWORK, LRWORK, M, N * .. * .. Array Arguments .. - DOUBLE COMPLEX A( LDA, * ), U( LDU, * ), V( LDV, * ), + COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), $ CWORK( LWORK ) DOUBLE PRECISION SVA( N ), RWORK( * ) INTEGER IWORK( * ) @@ -539,11 +539,11 @@ * .. Local Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) - DOUBLE COMPLEX CZERO, CONE + COMPLEX*16 CZERO, CONE PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ), CONE = ( 1.0D0, 0.0D0 ) ) * .. * .. Local Scalars .. - DOUBLE COMPLEX CTEMP + COMPLEX*16 CTEMP DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, $ COND_OK, CONDR1, CONDR2, ENTRA, ENTRAT, EPSLN, $ MAXPRJ, SCALEM, SCONDA, SFMIN, SMALL, TEMP1, @@ -554,18 +554,18 @@ $ NOSCAL, ROWPIV, RSVEC, TRANSP * .. * .. Intrinsic Functions .. - INTRINSIC ABS, DCMPLX, DCONJG, DLOG, DMAX1, DMIN1, DFLOAT, + INTRINSIC ABS, DCMPLX, DCONJG, DLOG, DMAX1, DMIN1, DBLE, $ MAX0, MIN0, NINT, DSQRT * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DZNRM2 - INTEGER IDAMAX + INTEGER IDAMAX, IZAMAX LOGICAL LSAME - EXTERNAL IDAMAX, LSAME, DLAMCH, DZNRM2 + EXTERNAL IDAMAX, IZAMAX, LSAME, DLAMCH, DZNRM2 * .. * .. External Subroutines .. - EXTERNAL ZCOPY, ZGELQF, ZGEQP3, ZGEQRF, ZLACPY, ZLASCL, - $ ZLASET, ZLASSQ, ZLASWP, ZUNGQR, ZUNMLQ, + EXTERNAL DLASSQ, ZCOPY, ZGELQF, ZGEQP3, ZGEQRF, ZLACPY, ZLASCL, + $ DLASCL, ZLASET, ZLASSQ, ZLASWP, ZUNGQR, ZUNMLQ, $ ZUNMQR, ZPOCON, DSCAL, ZDSCAL, ZSWAP, ZTRSM, XERBLA * EXTERNAL ZGESVJ @@ -640,7 +640,11 @@ * * Quick return for void matrix (Y3K safe) * #:) - IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN + IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN + IWORK(1:3) = 0 + RWORK(1:7) = 0 + RETURN + ENDIF * * Determine whether the matrix U should be M x N or M x M * @@ -665,7 +669,7 @@ * overflow. It is possible that this scaling pushes the smallest * column norm left from the underflow threshold (extreme case). * - SCALEM = ONE / DSQRT(DFLOAT(M)*DFLOAT(N)) + SCALEM = ONE / DSQRT(DBLE(M)*DBLE(N)) NOSCAL = .TRUE. GOSCAL = .TRUE. DO 1874 p = 1, N @@ -807,7 +811,7 @@ 1950 CONTINUE ELSE DO 1904 p = 1, M - RWORK(M+N+p) = SCALEM*ABS( A(p,IDAMAX(N,A(p,1),LDA)) ) + RWORK(M+N+p) = SCALEM*ABS( A(p,IZAMAX(N,A(p,1),LDA)) ) AATMAX = DMAX1( AATMAX, RWORK(M+N+p) ) AATMIN = DMIN1( AATMIN, RWORK(M+N+p) ) 1904 CONTINUE @@ -828,7 +832,7 @@ * XSC = ZERO TEMP1 = ONE - CALL ZLASSQ( N, SVA, 1, XSC, TEMP1 ) + CALL DLASSQ( N, SVA, 1, XSC, TEMP1 ) TEMP1 = ONE / TEMP1 * ENTRA = ZERO @@ -836,7 +840,7 @@ BIG1 = ( ( SVA(p) / XSC )**2 ) * TEMP1 IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1) 1113 CONTINUE - ENTRA = - ENTRA / DLOG(DFLOAT(N)) + ENTRA = - ENTRA / DLOG(DBLE(N)) * * Now, SVA().^2/Trace(A^* * A) is a point in the probability simplex. * It is derived from the diagonal of A^* * A. Do the same with the @@ -849,7 +853,7 @@ BIG1 = ( ( RWORK(p) / XSC )**2 ) * TEMP1 IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1) 1114 CONTINUE - ENTRAT = - ENTRAT / DLOG(DFLOAT(M)) + ENTRAT = - ENTRAT / DLOG(DBLE(M)) * * Analyze the entropies and decide A or A^*. Smaller entropy * usually means better input for the algorithm. @@ -905,9 +909,9 @@ * one should use ZGESVJ instead of ZGEJSV. * BIG1 = DSQRT( BIG ) - TEMP1 = DSQRT( BIG / DFLOAT(N) ) + TEMP1 = DSQRT( BIG / DBLE(N) ) * - CALL ZLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR ) + CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR ) IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN AAQQ = ( AAQQ / AAPP ) * TEMP1 ELSE @@ -1009,7 +1013,7 @@ * sigma_i < N*EPSLN*||A|| are flushed to zero. This is an * agressive enforcement of lower numerical rank by introducing a * backward error of the order of N*EPSLN*||A||. - TEMP1 = DSQRT(DFLOAT(N))*EPSLN + TEMP1 = DSQRT(DBLE(N))*EPSLN DO 3001 p = 2, N IF ( ABS(A(p,p)) .GE. (TEMP1*ABS(A(1,1))) ) THEN NR = NR + 1 @@ -1056,7 +1060,7 @@ TEMP1 = ABS(A(p,p)) / SVA(IWORK(p)) MAXPRJ = DMIN1( MAXPRJ, TEMP1 ) 3051 CONTINUE - IF ( MAXPRJ**2 .GE. ONE - DFLOAT(N)*EPSLN ) ALMORT = .TRUE. + IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE. END IF * * @@ -1136,7 +1140,7 @@ * IF ( L2PERT ) THEN * XSC = SQRT(SMALL) - XSC = EPSLN / DFLOAT(N) + XSC = EPSLN / DBLE(N) DO 4947 q = 1, NR CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO) DO 4949 p = 1, N @@ -1168,7 +1172,7 @@ * to drown denormals IF ( L2PERT ) THEN * XSC = SQRT(SMALL) - XSC = EPSLN / DFLOAT(N) + XSC = EPSLN / DBLE(N) DO 1947 q = 1, NR CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO) DO 1949 p = 1, NR @@ -1226,7 +1230,7 @@ CALL ZCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 ) CALL ZLACGV( NR-p+1, V(p,p), 1 ) 8998 CONTINUE - CALL ZLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) + CALL ZLASET('Upper', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV) * CALL ZGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U, $ LDU, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO ) @@ -1363,10 +1367,10 @@ CONDR1 = ONE / DSQRT(TEMP1) * .. here need a second oppinion on the condition number * .. then assume worst case scenario -* R1 is OK for inverse <=> CONDR1 .LT. DFLOAT(N) -* more conservative <=> CONDR1 .LT. SQRT(DFLOAT(N)) +* R1 is OK for inverse <=> CONDR1 .LT. DBLE(N) +* more conservative <=> CONDR1 .LT. SQRT(DBLE(N)) * - COND_OK = DSQRT(DSQRT(DFLOAT(NR))) + COND_OK = DSQRT(DSQRT(DBLE(NR))) *[TP] COND_OK is a tuning parameter. * IF ( CONDR1 .LT. COND_OK ) THEN @@ -1521,9 +1525,9 @@ CALL ZTRSM('L','U','C','N',NR,NR,CONE,CWORK(2*N+1), $ N,V,LDV) IF ( NR .LT. N ) THEN - CALL ZLASET('A',N-NR,NR,ZERO,CZERO,V(NR+1,1),LDV) - CALL ZLASET('A',NR,N-NR,ZERO,CZERO,V(1,NR+1),LDV) - CALL ZLASET('A',N-NR,N-NR,ZERO,CONE,V(NR+1,NR+1),LDV) + CALL ZLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV) + CALL ZLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV) + CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV) END IF CALL ZUNMQR('L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1), $ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR) @@ -1605,7 +1609,7 @@ * first QRF. Also, scale the columns to make them unit in * Euclidean norm. This applies to all cases. * - TEMP1 = DSQRT(DFLOAT(N)) * EPSLN + TEMP1 = DSQRT(DBLE(N)) * EPSLN DO 1972 q = 1, N DO 972 p = 1, N CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q) @@ -1635,7 +1639,7 @@ $ LDU, CWORK(N+1), LWORK-N, IERR ) * The columns of U are normalized. The cost is O(M*N) flops. - TEMP1 = DSQRT(DFLOAT(M)) * EPSLN + TEMP1 = DSQRT(DBLE(M)) * EPSLN DO 1973 p = 1, NR XSC = ONE / DZNRM2( M, U(1,p), 1 ) IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) @@ -1684,7 +1688,7 @@ DO 6972 p = 1, N CALL ZCOPY( N, CWORK(N+p), N, V(IWORK(p),1), LDV ) 6972 CONTINUE - TEMP1 = DSQRT(DFLOAT(N))*EPSLN + TEMP1 = DSQRT(DBLE(N))*EPSLN DO 6971 p = 1, N XSC = ONE / DZNRM2( N, V(1,p), 1 ) IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) @@ -1702,7 +1706,7 @@ END IF CALL ZUNMQR( 'Left', 'No Tr', M, N1, N, A, LDA, CWORK, U, $ LDU, CWORK(N+1), LWORK-N, IERR ) - TEMP1 = DSQRT(DFLOAT(M))*EPSLN + TEMP1 = DSQRT(DBLE(M))*EPSLN DO 6973 p = 1, N1 XSC = ONE / DZNRM2( M, U(1,p), 1 ) IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) @@ -1779,9 +1783,9 @@ NUMRANK = NINT(RWORK(2)) IF ( NR .LT. N ) THEN - CALL ZLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV ) - CALL ZLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV ) - CALL ZLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV ) + CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV ) + CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV ) + CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV ) END IF CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1), @@ -1791,7 +1795,7 @@ * first QRF. Also, scale the columns to make them unit in * Euclidean norm. This applies to all cases. * - TEMP1 = DSQRT(DFLOAT(N)) * EPSLN + TEMP1 = DSQRT(DBLE(N)) * EPSLN DO 7972 q = 1, N DO 8972 p = 1, N CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q) @@ -1836,7 +1840,7 @@ * Undo scaling, if necessary (and possible) * IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN - CALL ZLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR ) + CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR ) USCAL1 = ONE USCAL2 = ONE END IF