SUBROUTINE DSCLQ( M, N, A, LDA, L, LDL, WORK, IWORK, INFO ) IMPLICIT NONE INTEGER M, N, LDA, LDL, IWORK( * ) DOUBLE PRECISION A( LDA, * ), L( LDL, M ), WORK( M, * ) 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 CALL DRLQCQR( COUNT, M, N, A, LDA, WORK(1,1), WORK(1,M+1), $ WORK(1,2*M+1), IWORK, FLAG2, INFO ) IF (INFO .EQ. -1) THEN INFO = COUNT RETURN ENDIF IF (COUNT .EQ. 1) THEN !$OMP PARALLEL DO PRIVATE(J) DO I = 1, M DO J = 1, I-1 L(J, I) = ZERO ENDDO CALL DCOPY(M-(I-1),WORK(I,M+I),1,L(I,I),1) ENDDO !$OMP END PARALLEL DO ELSE CALL DTRMM('R','L','N','N',M,M,ONE,WORK(1,M+1),M,L,LDL) ENDIF IF (FLAG2 .EQV. .FALSE.) GO TO 10 RETURN END SUBROUTINE DSCLQ * SUBROUTINE DRLQCQR( COUNT, M, N, Q, LDQ, WORK, WORK2, WORK3, $ IWORK, FLAG2, INFO ) * IMPLICIT NONE * * .. Scalar Arguments .. INTEGER COUNT, M, N, LDQ, INFO, IWORK( * ) LOGICAL FLAG2 * .. * .. Array Arguments .. DOUBLE PRECISION Q( LDQ, * ), WORK( M, M ), WORK2( M, M ), $ 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 .. * INFO = 0 EPS = DLAMCH( 'Precision' ) OF = DLAMCH( 'O' ) I0 = 0 * FLAG2 = .TRUE. IF (COUNT .GE. 2) GO TO 10 !$OMP PARALLEL DO PRIVATE(TMP) DO K = 1, M TMP = DDOT( N, Q( K, 1 ), LDQ, Q( K, 1 ), LDQ) IF (TMP .GT. OF .OR. TMP .EQ. ZERO) THEN FLAG2 = .FALSE. ENDIF END DO !$OMP END PARALLEL DO * 10 IF (FLAG2) THEN CALL DSYRK('L','N',M,N,ONE,Q,LDQ,ZERO,WORK,M) * DO I = 1, M WORK3(I) = SQRT(WORK(I,I)) WORK3(M+I) = ONE/WORK3(I) ENDDO !$OMP PARALLEL DO PRIVATE(J) DO I = 1, M DO J = I, M WORK(J, I) = WORK(J, I)*WORK3(M+I)*WORK3(M+J) END DO END DO !$OMP END PARALLEL DO * SHIFT = ZERO * DO K = 1, 2*M !$OMP PARALLEL DO DO I = 1, M CALL DCOPY(M+1-I,WORK(I,I),1,WORK2(I,I),1) ENDDO !$OMP END PARALLEL DO IF (SHIFT .NE. ZERO) THEN DO I = 1, M WORK2(I,I) = WORK2(I,I) + SHIFT ENDDO ENDIF IINFO=0 CALL DPOTRF('L',M,WORK2,M,IINFO) IF (IINFO .EQ. 0) THEN IF (SHIFT .EQ. ZERO) THEN CALL DTRCON('I','L','N',M,WORK2,M,TMP,WORK3(M+1), $ IWORK,IINFO) IF (SQRT(SQRT(SQRT(EPS))) .GT. TMP) THEN FLAG2 = .FALSE. ELSE FLAG2 = .TRUE. END IF ELSE FLAG2 = .FALSE. END IF * !$OMP PARALLEL DO PRIVATE(J) DO I = 1, M DO J = I, M WORK2( J, I ) = WORK2( J, I )*WORK3(J) END DO END DO !$OMP END PARALLEL DO CALL DTRSM('L','L','N','N',M,N,ONE,WORK2,M,Q,LDQ) RETURN * ENDIF 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 (*,*) 'O',I0,SHIFT ENDDO * INFO = -1 RETURN * ELSE !$OMP PARALLEL DO DO K = 1, M WORK3(K) = DNRM2( N, Q( K, 1 ), LDQ ) CALL DSCAL( N, ONE/WORK3(K), Q( K, 1 ), LDQ ) END DO !$OMP END PARALLEL DO * CALL DSYRK('L','N',M,N,ONE,Q,LDQ,ZERO,WORK,M) * SHIFT = ZERO * DO K = 1, 2*M !$OMP PARALLEL DO DO I = 1, M CALL DCOPY(M+1-I,WORK(I,I),1,WORK2(I,I),1) ENDDO !$OMP END PARALLEL DO IF (SHIFT .NE. ZERO) THEN DO I = 1, M WORK2(I,I) = WORK2(I,I) + SHIFT ENDDO ENDIF IINFO=0 CALL DPOTRF('L',M,WORK2,M,IINFO) IF (IINFO .EQ. 0) THEN IF (SHIFT .EQ. ZERO) THEN CALL DTRCON('I','L','N',M,WORK2,M,TMP,WORK3(M+1), $ IWORK,IINFO) IF (SQRT(SQRT(SQRT(EPS))) .GT. TMP) THEN FLAG2 = .FALSE. ELSE FLAG2 = .TRUE. END IF ELSE FLAG2 = .FALSE. END IF * CALL DTRSM('L','L','N','N',M,N,ONE,WORK2,M,Q,LDQ) !$OMP PARALLEL DO PRIVATE(J) DO I = 1, M DO J = I, M WORK2( J, I ) = WORK2( J, I )*WORK3(J) END DO END DO !$OMP END PARALLEL DO RETURN * ENDIF 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 (*,*) 'O',I0,SHIFT ENDDO * INFO = -1 RETURN * ENDIF * END SUBROUTINE DRLQCQR