SUBROUTINE DSCQR( M, N, A, LDA, R, LDR, WORK, IWORK, INFO ) IMPLICIT NONE INTEGER M, N, LDA, LDR, IWORK( * ) DOUBLE PRECISION A( LDA, * ), R( LDR, N ), WORK( N, * ) DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE=1.0D+0, ZERO=0.0D+0 ) LOGICAL FLAG2 INTEGER COUNT, INFO, I, J COUNT = 0 10 COUNT = COUNT + 1 WRITE (*,*) COUNT CALL DRQRCQR( M, N, A, LDA, WORK(1,1), WORK(1,N+1), $ WORK(1,2*N+1), IWORK, FLAG2, INFO ) IF (INFO .EQ. -1) THEN INFO = COUNT RETURN ENDIF IF (COUNT .EQ. 1) THEN DO I = 1, N CALL DCOPY(I,WORK(1,N+I),1,R(1,I),1) DO J = I+1, N R(J, I) = ZERO ENDDO ENDDO ELSE CALL DTRMM('L','U','N','N',N,N,ONE,WORK(1,N+1),N,R,LDR) ENDIF IF (FLAG2 .EQV. .FALSE.) GO TO 10 RETURN END SUBROUTINE DSCQR * SUBROUTINE DSCQRH( COUNT, M, N, A, LDA, R, LDR, WORK, IWORK, $ INFO ) IMPLICIT NONE INTEGER COUNT, M, N, LDA, LDR, IWORK( * ) DOUBLE PRECISION A( LDA, * ), R( 5, LDR, N ), WORK( N, * ) DOUBLE PRECISION ZERO PARAMETER ( ZERO=0.0D+0 ) LOGICAL FLAG2 INTEGER INFO, I, J COUNT = 0 10 COUNT = COUNT + 1 IF (COUNT .EQ. 6) THEN INFO = COUNT RETURN ENDIF CALL DRQRCQR( M, N, A, LDA, WORK(1,1), WORK(1,N+1), $ WORK(1,2*N+1), IWORK, FLAG2, INFO ) IF (INFO .EQ. -1) RETURN DO I = 1, N CALL DCOPY(I,WORK(1,N+I),1,R(COUNT,1,I),1) DO J = I+1, N R(COUNT, J, I) = ZERO ENDDO ENDDO IF (FLAG2 .EQV. .FALSE.) GO TO 10 RETURN END SUBROUTINE DSCQRH * SUBROUTINE DRQRCQR( M, N, Q, LDQ, WORK, WORK2, WORK3, $ IWORK, FLAG2, INFO ) * IMPLICIT NONE * * .. Scalar Arguments .. INTEGER M, N, LDQ, IWORK( * ), INFO LOGICAL FLAG2 * .. * .. Array Arguments .. DOUBLE PRECISION Q( LDQ, * ), WORK( N, N ), WORK2( N, N ), $ WORK3( * ) * * ===================================================================== * * .. parameters .. DOUBLE PRECISION TWO, ONE, ZERO PARAMETER ( TWO=2.0D+0, ONE=1.0D+0, ZERO=0.0D+0 ) * .. * .. local scalars .. INTEGER IINFO, K, I, J, I0 DOUBLE PRECISION SHIFT, EPS, DLAMCH, TMP, OF, DNRM2, DDOT * .. * .. external functions .. EXTERNAL DSYRK, DPOTRF, DLAMCH, DNRM2, DDOT * .. * .. Executable Statements .. * EPS = DLAMCH( 'Precision' ) OF = DLAMCH( 'O' ) * CALL DSYRK('U','T',N,M,ONE,Q,LDQ,ZERO,WORK,N) DO K = 1, N IF (WORK(K,K) .GT. OF .OR. WORK(K,K) .EQ. ZERO) THEN GO TO 10 ENDIF END DO * SHIFT = ZERO I0 = 0 * DO K = 1, 2*N DO I = 1, N CALL DCOPY(I,WORK(1,I),1,WORK2(1,I),1) ENDDO IF (SHIFT .NE. ZERO) THEN DO I = 1, N WORK2(I,I) = WORK2(I,I) + SHIFT ENDDO ENDIF CALL DPOTRF('U',N,WORK2,N,IINFO) IF (IINFO .EQ. 0) THEN IF (SHIFT .EQ. ZERO) THEN CALL DTRCON('I','U','N',N,WORK2,N,TMP,WORK3(N+1), $ IWORK,IINFO) WRITE (*,*) 'Z',ONE/TMP IF (SQRT(SQRT(SQRT(EPS))) .GT. TMP) THEN FLAG2 = .FALSE. ELSE FLAG2 = .TRUE. END IF ELSE FLAG2 = .FALSE. END IF * CALL DTRSM('R','U','N','N',M,N,ONE,WORK2,N,Q,LDQ) INFO = 0 RETURN * ELSE IF (I0 .EQ. IINFO) THEN SHIFT = TWO*SHIFT ELSE I0 = IINFO TMP=SHIFT-WORK2(IINFO,IINFO) IF (TMP .NE. SHIFT) THEN SHIFT=TMP ELSE IF (SHIFT .EQ. ZERO) THEN TMP = ZERO DO I=1,N TMP = MAX(TMP,WORK(I,I)) ENDDO SHIFT=EPS*TMP ELSE SHIFT=(ONE+EPS)*SHIFT ENDIF ENDIF ENDIF WRITE (*,*) 'S',I0,SHIFT ENDIF ENDDO * INFO = -1 RETURN 10 DO K = 1, N WORK3(K) = DNRM2( M, Q( 1, K ), 1 ) CALL DSCAL( M, ONE/WORK3(K), Q( 1, K ), 1 ) END DO CALL DSYRK('U','T',N,M,ONE,Q,LDQ,ZERO,WORK,N) * SHIFT = ZERO I0 = 0 * DO K = 1, 2*N DO I = 1, N CALL DCOPY(I,WORK(1,I),1,WORK2(1,I),1) ENDDO IF (SHIFT .NE. ZERO) THEN DO I = 1, N WORK2(I,I) = WORK2(I,I) + SHIFT ENDDO ENDIF CALL DPOTRF('U',N,WORK2,N,IINFO) IF (IINFO .EQ. 0) THEN IF (SHIFT .EQ. ZERO) THEN CALL DTRCON('I','U','N',N,WORK2,N,TMP,WORK3(N+1), $ IWORK,IINFO) WRITE (*,*) 'Z',ONE/TMP IF (SQRT(SQRT(SQRT(EPS))) .GT. TMP) THEN FLAG2 = .FALSE. ELSE FLAG2 = .TRUE. END IF ELSE FLAG2 = .FALSE. END IF * CALL DTRSM('R','U','N','N',M,N,ONE,WORK2,N,Q,LDQ) DO I = 1, N CALL DSCAL( I, WORK3(I), WORK2( 1, I ), 1 ) END DO INFO = 0 RETURN * ELSE IF (I0 .EQ. IINFO) THEN SHIFT = TWO*SHIFT ELSE I0 = IINFO TMP=SHIFT-WORK2(IINFO,IINFO) IF (TMP .NE. SHIFT) THEN SHIFT=TMP ELSE IF (SHIFT .EQ. ZERO) THEN SHIFT=EPS ELSE SHIFT=(ONE+EPS)*SHIFT ENDIF ENDIF ENDIF WRITE (*,*) 'S',I0,SHIFT ENDIF ENDDO * INFO = -1 RETURN * END SUBROUTINE DRQRCQR