#! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'README' <<'END_OF_FILE' X *************************************************************************** X * All the software contained in this library is protected by copyright. * X * Permission to use, copy, modify, and distribute this software for any * X * purpose without fee is hereby granted, provided that this entire notice * X * is included in all copies of any software which is or includes a copy * X * or modification of this software and in all copies of the supporting * X * documentation for such software. * X *************************************************************************** X * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED * X * WARRANTY. IN NO EVENT, NEITHER THE AUTHORS, NOR THE PUBLISHER, NOR ANY * X * MEMBER OF THE EDITORIAL BOARD OF THE JOURNAL "NUMERICAL ALGORITHMS", * X * NOR ITS EDITOR-IN-CHIEF, BE LIABLE FOR ANY ERROR IN THE SOFTWARE, ANY * X * MISUSE OF IT OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF * X * USING THE SOFTWARE LIES WITH THE PARTY DOING SO. * X *************************************************************************** X * ANY USE OF THE SOFTWARE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE * X * ABOVE STATEMENT. * X *************************************************************************** X X AUTHOR: X X DARIO ANDREA BINI AND BEATRICE MEINI X UNIVERSITY OF PISA, ITALY X E-MAIL: bini@dm.unipi.it meini@dm.unipi.it X X REFERENCE: X X - IMPROVED CYCLIC REDUCTION FOR SOLVING QUEUING PROBLEMS X NUMERICAL ALGORITHMS, 15 (1997), PP. 57-74 X X SOFTWARE REVISION DATE: X X JANUARY 30, 1997 X X SOFTWARE LANGUAGE: X X FORTRAN 90 X X!=========================================================================! X! file README ! X!=========================================================================! X! IMPROVED CYCLIC REDUCTION FOR SOLVING QUEUEING PROBLEMS ! X! by D.A. Bini and B. Meini ! X! (bini@dm.unipi.it meini@dm.unipi.it) ! X! Fortran 90 Program version 1.0, January 30 1997 ! X! README FILE ! X!=========================================================================! X XThis package contains Fortran 90 programs for the numerical solution Xof the matrix equation: X X A_1+X A_2+X^2 A_3+....X^(m-1) A_m-X=O (1) X Xwhere X is the (nb x nb) unknown matrix and the (nb x nb) nonnegative Xmatrices A_i, i=1,2,..., are such that A_1+A_2+...+A_m is column Xstochastic. X X X XThe package contains the following files: X XREADME : this file Xpwcr_drv.f90 : driver program Xpwcr_sub.f90 : main subroutines Xpwcr_int.f90 : interfaces Xpwcr_fft.f90 : fft subroutines Xfft_int.f90 : fft interfaces Xkgesv1.f90 : Lapack fortran 90 subroutines Xlaauxmod.f90 : Lapack interfaces Xlaf90mod.f90 : Lapack interfaces Xlaf77mod.f90 : Lapack interfaces Xsolve.f : Lapack fortran 77 subroutines XMakefile : make file Xdata : input data file Xresults : output file generated by the driver X X-------------------------------------------------------------------------- X XIn order to create the executable file pwcr (Point-Wise Cyclic Reduction), Xit is sufficient to type `make' X XDuring the compilation of the file Xsolve.f, you could receive some warning messages, like: X XWarning: solve.f, line 1324: Unused dummy argument N3 X detected at END@ XWarning: solve.f, line 1324: Unused dummy argument OPTS X detected at END@ XExtension: solve.f, line xxxx: Byte count on numeric data type X detected at *@16 X XThese warning messages, originated from the Lapack fortran 77 subroutines, Xdo not affect the correctness of the compilation. X XOnce compilation has been performed, the user may remove the files *.o Xand *.mod. X XIn order to check the correct installation of the package Xrun the executable `pwcr'. XAt the request `input file name' type `data'. XIn this way the file data.out will be created. This file, up to within Xsmall numerical differences due to a possibly different floating point Xarithmetic, should coincide with the file `results'. X X------------------------------------------------------------------ X XIn order to use the program `pwcr' with different input file, the Xuser should create her/his own file, say, `my_file' organized Xin the following way: X XEach row must contain the following data: X X1-st row: the dimension of the blocks ( integer ) X2-nd row: the number of blocks A_i ( integer ) X3-rd row: the precision level ( real(kind(0.d0)) ), say, 1.d0d-12 XThe subsequent rows: the values L, I, J, VAL, where X L : (integer) is the index of the block A_L X I,J : (integers) are the indices of the (I,J)-th X entry of the block A_L X VAL : ( real(kind(0.d0)) ) is the value of the X (I,J)-th entry of the block A_L XThe last row: the quartuple 0,0,0,0.0d0 that denotes X the end of the file X XFor the values L, I, J, VAL, that are not reported in this file it is Xassumed VAL=0.d0. X X XThe file `my_file.out' created by the program contains the output Xorganized in the following way: X X --The residual error ||A_1+G A_2+G^2 A_3+....G^(m-1) A_m-G ||_1, where X G is the computed approximation of the solution of (1) and X ||.||_1 denotes the 1-norm. X --The entries of the matrix G arranged row-wise X X----------------------------------------------------------------------- X END_OF_FILE if test 5932 -ne `wc -c <'README'`; then echo shar: \"'README'\" unpacked with wrong size! fi # end of 'README' fi if test -f 'Makefile' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'Makefile'\" else echo shar: Extracting \"'Makefile'\" \(1448 characters\) sed "s/^X//" >'Makefile' <<'END_OF_FILE' X#=========================================================================# X# file Makefile # X#=========================================================================# X# IMPROVED CYCLIC REDUCTION FOR SOLVING QUEUEING PROBLEMS # X# by D.A. Bini and B. Meini # X# (bini@dm.unipi.it meini@dm.unipi.it) # X# Fortran 90 Program version 1.0, January 30, 1997 # X# Makefile # X#=========================================================================# X XF90 = f90 XFFLAGS = -O X XOBJ = pwcr_int.o fft_int.o laauxmod.o laf77mod.o laf90mod.o pwcr_sub.o pwcr_fft.o kgesv1.o solve.o pwcr_drv.o X Xall: pwcr X Xpwcr:$(OBJ) X $(F90) $(OBJ) -o pwcr X Xpwcr_int.o:pwcr_int.f90 X $(F90) $(FFLAGS) -c pwcr_int.f90 X Xfft_int.o:fft_int.f90 X $(F90) $(FFLAGS) -c fft_int.f90 X Xlaf77mod.o:laf77mod.f90 X $(F90) $(FFLAGS) -c laf77mod.f90 X Xlaf90mod.o:laf90mod.f90 X $(F90) $(FFLAGS) -c laf90mod.f90 X Xlaauxmod.o:laauxmod.f90 X $(F90) $(FFLAGS) -c laauxmod.f90 X Xpwcr_drv.o:pwcr_drv.f90 X $(F90) $(FFLAGS) -c pwcr_drv.f90 X Xpwcr_sub.o:pwcr_sub.f90 X $(F90) $(FFLAGS) -c pwcr_sub.f90 X Xpwcr_fft.o:pwcr_fft.f90 X $(F90) $(FFLAGS) -c pwcr_fft.f90 X Xkgesv1.o:kgesv1.f90 X $(F90) $(FFLAGS) -c kgesv1.f90 X Xsolve.o:solve.f X $(F90) $(FFLAGS) -c solve.f X X X X X X X X X X END_OF_FILE if test 1448 -ne `wc -c <'Makefile'`; then echo shar: \"'Makefile'\" unpacked with wrong size! fi # end of 'Makefile' fi if test -f 'fft_int.f90' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'fft_int.f90'\" else echo shar: Extracting \"'fft_int.f90'\" \(5452 characters\) sed "s/^X//" >'fft_int.f90' <<'END_OF_FILE' X!=========================================================================! X! file fft_int.f90 ! X!=========================================================================! X! IMPROVED CYCLIC REDUCTION FOR SOLVING QUEUEING PROBLEMS ! X! by D.A. Bini and B. Meini ! X! (bini@dm.unipi.it meini@dm.unipi.it) ! X! Fortran 90 Program version 1.0, January 30 1997 ! X! Interface File for FFT subroutines ! X!=========================================================================! X! Interface file for the following FFT subroutines: ! X! iffts1 ! X! ffts2 ! X! iffts2 ! X! ffts1 ! X! fillroots ! X! fft1 ! X! ifft1 ! X! ftb1 ! X! ftb2 ! X! iftb2 ! X! iftb1 ! X! twiddle ! X! itwiddle ! X!=========================================================================! XMODULE fft_interface X INTERFACE X SUBROUTINE iffts1(u,v) X IMPLICIT NONE X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: u, v X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: wr, wi X COMMON wr, wi X END SUBROUTINE iffts1 X END INTERFACE X !================================================================= X X INTERFACE X SUBROUTINE ffts2(u,v) X IMPLICIT NONE X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: u, v X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: wr, wi X COMMON wr, wi X END SUBROUTINE ffts2 X END INTERFACE X !================================================================= X X INTERFACE X SUBROUTINE iffts2(u,v,t) X IMPLICIT NONE X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: u, v,t X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: wr, wi X COMMON wr, wi X END SUBROUTINE iffts2 X END INTERFACE X !================================================================= X X INTERFACE X SUBROUTINE ffts1(u,v) X IMPLICIT NONE X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: u, v X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: wr, wi X COMMON wr, wi X END SUBROUTINE ffts1 X END INTERFACE X !================================================================= X X INTERFACE X SUBROUTINE fillroots(n) X IMPLICIT NONE X REAL(KIND(0.d0)),DIMENSION(:),ALLOCATABLE,SAVE :: wwr, wwi X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: wr, wi X INTEGER :: n X COMMON wr, wi X END SUBROUTINE fillroots X END INTERFACE X !================================================================= X X INTERFACE X SUBROUTINE fft1(x,y) X IMPLICIT NONE X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: x,y X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: wr, wi X COMMON wr, wi X END SUBROUTINE FFT1 X END INTERFACE X !================================================================= X X INTERFACE X SUBROUTINE ifft1(x,y) X IMPLICIT NONE X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: x,y X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: wr, wi X COMMON wr, wi X END SUBROUTINE ifft1 X END INTERFACE X !================================================= X X INTERFACE X SUBROUTINE ftb1(a,ta) X IMPLICIT NONE X REAL(KIND(0.d0)),DIMENSION(:,:,:),POINTER :: a, ta X END SUBROUTINE ftb1 X END INTERFACE X !===================================================== X X INTERFACE X SUBROUTINE ftb2(a,ta) X IMPLICIT NONE X REAL(KIND(0.d0)),DIMENSION(:,:,:),POINTER :: a, ta X END SUBROUTINE ftb2 X END INTERFACE X !===================================================== X X INTERFACE X SUBROUTINE iftb2(ta,a) X IMPLICIT NONE X REAL(KIND(0.d0)),DIMENSION(:,:,:),POINTER :: a,ta X END SUBROUTINE iftb2 X END INTERFACE X !===================================================== X X INTERFACE X SUBROUTINE iftb1(ta,a) X IMPLICIT NONE X REAL(KIND(0.d0)),DIMENSION(:,:,:),POINTER :: a, ta X END SUBROUTINE iftb1 X END INTERFACE X !===================================================== X X INTERFACE X SUBROUTINE twiddle(zr,zi) X IMPLICIT NONE X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: zr,zi X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: wr, wi X COMMON wr, wi X END SUBROUTINE twiddle X END INTERFACE X !=================================================== X X INTERFACE X SUBROUTINE itwiddle(zr,zi) X IMPLICIT NONE X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: zr, zi X REAL(KIND(0.d0)),DIMENSION(:),POINTER :: wr, wi X COMMON wr, wi X END SUBROUTINE itwiddle X END INTERFACE X XEND MODULE fft_interface X X X X X X X X X X X X END_OF_FILE if test 5452 -ne `wc -c <'fft_int.f90'`; then echo shar: \"'fft_int.f90'\" unpacked with wrong size! fi # end of 'fft_int.f90' fi if test -f 'kgesv1.f90' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'kgesv1.f90'\" else echo shar: Extracting \"'kgesv1.f90'\" \(6241 characters\) sed "s/^X//" >'kgesv1.f90' <<'END_OF_FILE' X!================================================================ X! file kgesv1.f90 X!================================================================ X! /netlib/lapack90 Preliminary version November, 1996 X! X! file: kgesv.f90 X! for: Fortran 90 wrapper for LAPACK routines xGESV X!================================================================ X X SUBROUTINE DGESV_f90(A,B,IPIV,INFO) X! .. Use Statements .. X USE LA_PRECISION, ONLY: WP => DP X USE LA_AUX, ONLY: ERINFO X USE LAPACK77_INTERFACES, ONLY: GESV_F77 => DGESV X! .. Implicit Statement .. X IMPLICIT NONE X! .. Scalar Arguments .. X INTEGER, INTENT(OUT), OPTIONAL :: INFO X! .. Array Arguments .. X INTEGER, INTENT(OUT), OPTIONAL, TARGET :: IPIV(:) X REAL(WP), INTENT(INOUT) :: A(:,:), B(:,:) X! .. Parameters .. X CHARACTER(LEN=7), PARAMETER :: SRNAME = 'LA_GESV' X! .. Local Scalars .. X INTEGER :: LD, LINFO, NRHS, N X! .. Local Pointers .. X INTEGER, POINTER :: LPIV(:) X! .. Intrinsic Functions .. X! INTRINSIC ALLOCATE, DEALLOCATE, MAX, PRESENT, SIZE X INTRINSIC MAX, PRESENT, SIZE X! X! .. Executable Statements .. X LINFO = 0 X N = SIZE(A, 1) X IF( SIZE( A, 2 ) /= N ) THEN X LINFO = -1 X ELSE IF( SIZE( B, 1 ) /= N ) THEN X LINFO = -2 X ELSE X IF( PRESENT( IPIV ) )THEN X IF( SIZE( IPIV ) /= N ) LINFO = -3 X END IF X END IF X! X IF ( LINFO == 0 ) THEN X LD = MAX( 1, N ) X NRHS = SIZE(B,2) X IF( PRESENT( IPIV ) )THEN X LPIV => IPIV X ELSE X ALLOCATE(LPIV(N)) X END IF X CALL GESV_F77( N, NRHS, A, LD, LPIV, B, LD, LINFO ) X IF(.NOT.PRESENT(IPIV)) DEALLOCATE(LPIV) X END IF X CALL ERINFO(LINFO,SRNAME,INFO) X! X END SUBROUTINE DGESV_F90 X X SUBROUTINE ZGESV_f90(A,B,IPIV,INFO) X! .. Use Statements .. X USE LA_PRECISION, ONLY: WP => DP X USE LA_AUX, ONLY: ERINFO X USE LAPACK77_INTERFACES, ONLY: GESV_F77 => ZGESV X! .. Implicit Statement .. X IMPLICIT NONE X! .. Scalar Arguments .. X INTEGER, INTENT(OUT), OPTIONAL :: INFO X! .. Array Arguments .. X INTEGER, INTENT(OUT), OPTIONAL, TARGET :: IPIV(:) X COMPLEX(WP), INTENT(INOUT) :: A(:,:), B(:,:) X! .. Parameters .. X CHARACTER(LEN=7), PARAMETER :: SRNAME = 'LA_GESV' X! .. Local Scalars .. X INTEGER :: LD, LINFO, NRHS, N X! .. Local Pointers .. X INTEGER, POINTER :: LPIV(:) X! .. Intrinsic Functions .. X! INTRINSIC ALLOCATE, DEALLOCATE, MAX, PRESENT, SIZE X INTRINSIC MAX, PRESENT, SIZE X! X! .. Executable Statements .. X LINFO = 0 X N = SIZE(A, 1) X IF( SIZE( A, 2 ) /= N ) THEN X LINFO = -1 X ELSE IF( SIZE( B, 1 ) /= N ) THEN X LINFO = -2 X ELSE X IF( PRESENT( IPIV ) )THEN X IF( SIZE( IPIV ) /= N ) LINFO = -3 X END IF X END IF X IF ( LINFO == 0 ) THEN X LD = MAX( 1, N ) X NRHS = SIZE(B,2) X IF( PRESENT( IPIV ) )THEN X LPIV => IPIV X ELSE X ALLOCATE(LPIV(N)) X END IF X CALL GESV_F77( N, NRHS, A, LD, LPIV, B, LD, LINFO ) X IF(.NOT.PRESENT(IPIV)) DEALLOCATE(LPIV) X END IF X CALL ERINFO(LINFO,SRNAME,INFO) X END SUBROUTINE ZGESV_F90 X X SUBROUTINE D1GESV_f90(A,B,IPIV,INFO) X! .. Use Statements .. X USE LA_PRECISION, ONLY: WP => DP X USE LA_AUX, ONLY: ERINFO X USE LAPACK77_INTERFACES, ONLY: GESV_F77 => DGESV X! .. Implicit Statement .. X IMPLICIT NONE X! .. Scalar Arguments .. X INTEGER, INTENT(OUT), OPTIONAL :: INFO X! .. Array Arguments .. X INTEGER, INTENT(OUT), OPTIONAL, TARGET :: IPIV(:) X REAL(WP), INTENT(INOUT) :: A(:,:), B(:) X! .. Parameters .. X CHARACTER(LEN=7), PARAMETER :: SRNAME = 'LA_GESV' X! .. Local Scalars .. X INTEGER :: LD, LINFO, NRHS, N X! .. Local Pointers .. X INTEGER, POINTER :: LPIV(:) X! .. Intrinsic Functions .. X! INTRINSIC ALLOCATE, DEALLOCATE, MAX, PRESENT, SIZE X INTRINSIC MAX, PRESENT, SIZE X! X! .. Executable Statements .. X LINFO = 0 X N = SIZE(A, 1) X IF( SIZE( A, 2 ) /= N ) THEN X LINFO = -1 X ELSE IF( SIZE( B ) /= N ) THEN X LINFO = -2 X ELSE X IF( PRESENT( IPIV ) )THEN X IF( SIZE( IPIV ) /= N ) LINFO = -3 X END IF X END IF X IF ( LINFO == 0 ) THEN X LD = MAX( 1, N ) X NRHS = 1 X IF( PRESENT( IPIV ) )THEN X LPIV => IPIV X ELSE X ALLOCATE(LPIV(N)) X END IF X CALL GESV_F77( N, NRHS, A, LD, LPIV, B, LD, LINFO ) X IF(.NOT.PRESENT(IPIV)) DEALLOCATE(LPIV) X END IF X CALL ERINFO(LINFO,SRNAME,INFO) X END SUBROUTINE D1GESV_F90 X X SUBROUTINE Z1GESV_f90(A,B,IPIV,INFO) X! .. Use Statements .. X USE LA_PRECISION, ONLY: WP => DP X USE LA_AUX, ONLY: ERINFO X USE LAPACK77_INTERFACES, ONLY: GESV_F77 => ZGESV X! .. Implicit Statement .. X IMPLICIT NONE X! .. Scalar Arguments .. X INTEGER, INTENT(OUT), OPTIONAL :: INFO X! .. Array Arguments .. X INTEGER, INTENT(OUT), OPTIONAL, TARGET :: IPIV(:) X COMPLEX(WP), INTENT(INOUT) :: A(:,:), B(:) X! .. Parameters .. X CHARACTER(LEN=7), PARAMETER :: SRNAME = 'LA_GESV' X! .. Local Scalars .. X INTEGER :: LD, LINFO, NRHS, N X! .. Local Pointers .. X INTEGER, POINTER :: LPIV(:) X! .. Intrinsic Functions .. X! INTRINSIC ALLOCATE, DEALLOCATE, MAX, PRESENT, SIZE X INTRINSIC MAX, PRESENT, SIZE X! X! .. Executable Statements .. X LINFO = 0 X N = SIZE(A, 1) X IF( SIZE( A, 2 ) /= N ) THEN X LINFO = -1 X ELSE IF( SIZE( B ) /= N ) THEN X LINFO = -2 X ELSE X IF( PRESENT( IPIV ) )THEN X IF( SIZE( IPIV ) /= N ) LINFO = -3 X END IF X END IF X IF ( LINFO == 0 ) THEN X LD = MAX( 1, N ) X NRHS = 1 X IF( PRESENT( IPIV ) )THEN X LPIV => IPIV X ELSE X ALLOCATE(LPIV(N)) X END IF X CALL GESV_F77( N, NRHS, A, LD, LPIV, B, LD, LINFO ) X IF(.NOT.PRESENT(IPIV)) DEALLOCATE(LPIV) X END IF X CALL ERINFO(LINFO,SRNAME,INFO) X END SUBROUTINE Z1GESV_F90 END_OF_FILE if test 6241 -ne `wc -c <'kgesv1.f90'`; then echo shar: \"'kgesv1.f90'\" unpacked with wrong size! fi # end of 'kgesv1.f90' fi if test -f 'laauxmod.f90' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'laauxmod.f90'\" else echo shar: Extracting \"'laauxmod.f90'\" \(3407 characters\) sed "s/^X//" >'laauxmod.f90' <<'END_OF_FILE' X!================================================================ X! file laauxmod.f90 X!================================================================ X! /netlib/lapack90 Preliminary version November, 1996 X! X! file: laauxmod.f90 X! for: auxiliary routines X!================================================================ X X XMODULE LA_PRECISION X! X! Defines single and double precision parameters, sp and dp. X! These values are compiler dependent. X! X INTEGER, PARAMETER :: SP=KIND(1.0), DP=KIND(1.0D0) X! XEND MODULE LA_PRECISION X XMODULE LA_AUX X! XCONTAINS X! XLOGICAL FUNCTION LSAME( CA, CB ) X! X! -- LAPACK auxiliary routine (version 2.0) -- X! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., X! Courant Institute, Argonne National Lab, and Rice University X! September 30, 1994 X! X! Purpose X! ======= X! X! LSAME tests if CA is the same letter as CB regardless of case. X! X! Parameters X! ========== X! X! CA (input) CHARACTER*1 X! CB (input) CHARACTER*1 X! Characters to be compared. X! X! .. Scalar Arguments .. X CHARACTER*1, INTENT(IN) :: CA, CB X! .. Parameters .. X INTEGER, PARAMETER :: IOFF=32 X! .. Local Scalars .. X INTEGER :: INTA, INTB, ZCODE X! .. Intrinsic Functions .. X INTRINSIC ICHAR X! X! .. Executable Statements .. X! X! Test if the characters are equal X! X LSAME = CA == CB X! X! Now test for equivalence X! X IF ( .NOT.LSAME ) THEN X! X! Use 'Z' rather than 'A' so that ASCII can be detected on Prime X! machines, on which ICHAR returns a value with bit 8 set. X! ICHAR('A') on Prime machines returns 193 which is the same as X! ICHAR('A') on an EBCDIC machine. X! X ZCODE = ICHAR( 'Z' ) X! X INTA = ICHAR( CA ) X INTB = ICHAR( CB ) X! X IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN X! X! ASCII is assumed - ZCODE is the ASCII code of either lower or X! upper case 'Z'. X! X IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 X IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 X! X ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN X! X! EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or X! upper case 'Z'. X! X IF( INTA.GE.129 .AND. INTA.LE.137 .OR. & X! INTA.GE.145 .AND. INTA.LE.153 .OR. & X INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 X IF( INTB.GE.129 .AND. INTB.LE.137 .OR. & X INTB.GE.145 .AND. INTB.LE.153 .OR. & X INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 X! X ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN X! X! ASCII is assumed, on Prime machines - ZCODE is the ASCII code X! plus 128 of either lower or upper case 'Z'. X! X IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 X IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 X END IF X LSAME = INTA == INTB X END IF XEND FUNCTION LSAME X XSUBROUTINE ERINFO(LINFO, SRNAME, INFO) X! .. Scalar Arguments .. X CHARACTER( LEN = * ), INTENT(IN) :: SRNAME X INTEGER , INTENT(IN) :: LINFO X INTEGER , INTENT(INOUT), OPTIONAL :: INFO X! X! .. Executable Statements .. X! X IF( PRESENT(INFO) ) INFO = LINFO X IF( LINFO < 0 .OR. LINFO>0 .AND. .NOT.PRESENT(INFO) )THEN X WRITE (*,*) 'Program terminated in LAPACK_90 subroutine ', SRNAME X WRITE (*,*) 'Error indicator, INFO = ', LINFO X STOP X END IF XEND SUBROUTINE ERINFO X! XEND MODULE LA_AUX END_OF_FILE if test 3407 -ne `wc -c <'laauxmod.f90'`; then echo shar: \"'laauxmod.f90'\" unpacked with wrong size! fi # end of 'laauxmod.f90' fi if test -f 'laf77mod.f90' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'laf77mod.f90'\" else echo shar: Extracting \"'laf77mod.f90'\" \(16245 characters\) sed "s/^X//" >'laf77mod.f90' <<'END_OF_FILE' X!================================================================ X! file laf77mod.f90 X!================================================================ X! /netlib/lapack90 Preliminary version November, 1996 X! X! file: laf77mod.f90 X! for: Interface module for the LAPACK Fortran 77 routines X!================================================================ X XMODULE LAPACK77_INTERFACES X INTERFACE X! X SUBROUTINE SGESV( N, NRHS, A, LDA, PIV, B, LDB, INFO ) X USE LA_PRECISION, ONLY: SP X INTEGER, INTENT(IN) :: LDA, LDB, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(OUT) :: PIV(*) X REAL(SP), INTENT(INOUT) :: A(LDA,*), B(LDB,NRHS) X END SUBROUTINE SGESV X SUBROUTINE DGESV( N, NRHS, A, LDA, PIV, B, LDB, INFO ) X USE LA_PRECISION, ONLY: DP X INTEGER, INTENT(IN) :: LDA, LDB, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(OUT) :: PIV(*) X REAL(DP), INTENT(INOUT) :: A(LDA,*), B(LDB,NRHS) X END SUBROUTINE DGESV X SUBROUTINE CGESV( N, NRHS, A, LDA, PIV, B, LDB, INFO ) X USE LA_PRECISION, ONLY: SP X INTEGER, INTENT(IN) :: LDA, LDB, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(OUT) :: PIV(*) X COMPLEX(SP), INTENT(INOUT) :: A(LDA,*), B(LDB,NRHS) X END SUBROUTINE CGESV X SUBROUTINE ZGESV( N, NRHS, A, LDA, PIV, B, LDB, INFO ) X USE LA_PRECISION, ONLY: DP X INTEGER, INTENT(IN) :: LDA, LDB, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(OUT) :: PIV(*) X COMPLEX(DP), INTENT(INOUT) :: A(LDA,*), B(LDB,NRHS) X END SUBROUTINE ZGESV X! X SUBROUTINE SGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, & X PIV, EQUED, R, C, B, LDB, X, LDX, & X RCOND, FERR, BERR, WORK, IWORK, INFO ) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN) :: TRANS, FACT X CHARACTER(LEN=1), INTENT(INOUT) :: EQUED X INTEGER, INTENT(IN) :: LDA, LDAF, LDB, LDX, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(OUT) :: IWORK(*) X INTEGER, INTENT(INOUT) :: PIV(*) X REAL(WP), INTENT(OUT) :: RCOND X REAL(WP), INTENT(OUT) :: FERR(*), BERR(*) X REAL(WP), INTENT(OUT) :: X(LDX,*), WORK(*) X REAL(WP), INTENT(INOUT) :: R(*), C(*) X REAL(WP), INTENT(INOUT) :: A(LDA,*), AF(LDAF,*), B(LDB,*) X END SUBROUTINE SGESVX X SUBROUTINE DGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, & X PIV, EQUED, R, C, B, LDB, X, LDX, & X RCOND, FERR, BERR, WORK, IWORK, INFO ) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN) :: TRANS, FACT X CHARACTER(LEN=1), INTENT(INOUT) :: EQUED X INTEGER, INTENT(IN) :: LDA, LDAF, LDB, LDX, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(OUT) :: IWORK(*) X INTEGER, INTENT(INOUT) :: PIV(*) X REAL(WP), INTENT(OUT) :: RCOND X REAL(WP), INTENT(OUT) :: FERR(*), BERR(*) X REAL(WP), INTENT(OUT) :: X(LDX,*), WORK(*) X REAL(WP), INTENT(INOUT) :: R(*), C(*) X REAL(WP), INTENT(INOUT) :: A(LDA,*), AF(LDAF,*), B(LDB,*) X END SUBROUTINE DGESVX X SUBROUTINE CGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, & X PIV, EQUED, R, C, B, LDB, X, LDX, & X RCOND, FERR, BERR, WORK, RWORK, INFO ) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN) :: TRANS, FACT X CHARACTER(LEN=1), INTENT(INOUT) :: EQUED X INTEGER, INTENT(IN) :: LDA, LDAF, LDB, LDX, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(INOUT) :: PIV(*) X REAL(WP), INTENT(OUT) :: RCOND X REAL(WP), INTENT(OUT) :: FERR(*), BERR(*), RWORK(*) X REAL(WP), INTENT(INOUT) :: R(*), C(*) X COMPLEX(WP), INTENT(OUT) :: X(LDX,*), WORK(*) X COMPLEX(WP), INTENT(INOUT) :: A(LDA,*), AF(LDAF,*), B(LDB,*) X END SUBROUTINE CGESVX X SUBROUTINE ZGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, & X PIV, EQUED, R, C, B, LDB, X, LDX, & X RCOND, FERR, BERR, WORK, RWORK, INFO ) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN) :: TRANS, FACT X CHARACTER(LEN=1), INTENT(INOUT) :: EQUED X INTEGER, INTENT(IN) :: LDA, LDAF, LDB, LDX, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(INOUT) :: PIV(*) X REAL(WP), INTENT(OUT) :: RCOND X REAL(WP), INTENT(OUT) :: FERR(*), BERR(*), RWORK(*) X REAL(WP), INTENT(INOUT) :: R(*), C(*) X COMPLEX(WP), INTENT(OUT) :: X(LDX,*), WORK(*) X COMPLEX(WP), INTENT(INOUT) :: A(LDA,*), AF(LDAF,*), B(LDB,*) X END SUBROUTINE ZGESVX X! X SUBROUTINE SGETRF( M, N, A, LDA, PIV, INFO ) X USE LA_PRECISION, ONLY: SP X INTEGER, INTENT(IN) :: LDA, M, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT( OUT ) :: PIV( * ) X REAL(SP), INTENT( INOUT ) :: A( LDA, * ) X END SUBROUTINE SGETRF X SUBROUTINE DGETRF( M, N, A, LDA, PIV, INFO ) X USE LA_PRECISION, ONLY: DP X INTEGER, INTENT(IN) :: LDA, M, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT( OUT ) :: PIV( * ) X REAL(DP), INTENT( INOUT ) :: A( LDA, * ) X END SUBROUTINE DGETRF X SUBROUTINE CGETRF( M, N, A, LDA, PIV, INFO ) X USE LA_PRECISION, ONLY: SP X INTEGER, INTENT(IN) :: LDA, M, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT( OUT ) :: PIV( * ) X COMPLEX(SP), INTENT( INOUT ) :: A( LDA, * ) X END SUBROUTINE CGETRF X SUBROUTINE ZGETRF( M, N, A, LDA, PIV, INFO ) X USE LA_PRECISION, ONLY: DP X INTEGER, INTENT(IN) :: LDA, M, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT( OUT ) :: PIV( * ) X COMPLEX(DP), INTENT( INOUT ) :: A( LDA, * ) X END SUBROUTINE ZGETRF X! X SUBROUTINE SGETRS( TRANS, N, NRHS, A, LDA, PIV, B, LDB, INFO ) X USE LA_PRECISION, ONLY: SP X CHARACTER(LEN=1), INTENT(IN) :: TRANS X INTEGER, INTENT(IN) :: LDA, LDB, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(IN) :: PIV(*) X REAL(SP), INTENT(IN) :: A(LDA,*) X REAL(SP), INTENT(INOUT) :: B(LDB,NRHS) X END SUBROUTINE SGETRS X SUBROUTINE DGETRS( TRANS, N, NRHS, A, LDA, PIV, B, LDB, INFO ) X USE LA_PRECISION, ONLY: DP X CHARACTER(LEN=1), INTENT(IN) :: TRANS X INTEGER, INTENT(IN) :: LDA, LDB, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(IN) :: PIV(*) X REAL(DP), INTENT(IN) :: A(LDA,*) X REAL(DP), INTENT(INOUT) :: B(LDB,NRHS) X END SUBROUTINE DGETRS X SUBROUTINE CGETRS( TRANS, N, NRHS, A, LDA, PIV, B, LDB, INFO ) X USE LA_PRECISION, ONLY: SP X CHARACTER(LEN=1), INTENT(IN) :: TRANS X INTEGER, INTENT(IN) :: LDA, LDB, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(IN) :: PIV(*) X COMPLEX(SP), INTENT(IN) :: A(LDA,*) X COMPLEX(SP), INTENT(INOUT) :: B(LDB,NRHS) X END SUBROUTINE CGETRS X SUBROUTINE ZGETRS( TRANS, N, NRHS, A, LDA, PIV, B, LDB, INFO ) X USE LA_PRECISION, ONLY: DP X CHARACTER(LEN=1), INTENT(IN) :: TRANS X INTEGER, INTENT(IN) :: LDA, LDB, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(IN) :: PIV(*) X COMPLEX(DP), INTENT(IN) :: A(LDA,*) X COMPLEX(DP), INTENT(INOUT) :: B(LDB,NRHS) X END SUBROUTINE ZGETRS X! X SUBROUTINE SGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) X USE LA_PRECISION, ONLY: SP X INTEGER, INTENT(IN) :: LDA, LWORK, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(IN) :: IPIV(*) X REAL(SP), INTENT(OUT) :: WORK(*) X REAL(SP), INTENT(INOUT) :: A(LDA,*) X END SUBROUTINE SGETRI X SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) X USE LA_PRECISION, ONLY: DP X INTEGER, INTENT(IN) :: LDA, LWORK, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(IN) :: IPIV(*) X REAL(DP), INTENT(OUT) :: WORK(*) X REAL(DP), INTENT(INOUT) :: A(LDA,*) X END SUBROUTINE DGETRI X SUBROUTINE CGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) X USE LA_PRECISION, ONLY: SP X INTEGER, INTENT(IN) :: LDA, LWORK, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(IN) :: IPIV(*) X COMPLEX(SP), INTENT(OUT) :: WORK(*) X COMPLEX(SP), INTENT(INOUT) :: A(LDA,*) X END SUBROUTINE CGETRI X SUBROUTINE ZGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) X USE LA_PRECISION, ONLY: DP X INTEGER, INTENT(IN) :: LDA, LWORK, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(IN) :: IPIV(*) X COMPLEX(DP), INTENT(OUT) :: WORK(*) X COMPLEX(DP), INTENT(INOUT) :: A(LDA,*) X END SUBROUTINE ZGETRI X! X SUBROUTINE SGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, PIV, B, & X LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO ) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN) :: TRANS X INTEGER, INTENT(IN) :: LDA, LDAF, LDB, LDX, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(IN) :: PIV(*) X INTEGER, INTENT(OUT) :: IWORK(*) X REAL(WP), INTENT(OUT) :: FERR(*), BERR(*), WORK(*) X REAL(WP), INTENT(IN) :: A(LDA,*), AF(LDAF,*), B(LDB,*) X REAL(WP), INTENT(INOUT) :: X(LDX,*) X END SUBROUTINE SGERFS X SUBROUTINE DGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, PIV, B, & X LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO ) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN) :: TRANS X INTEGER, INTENT(IN) :: LDA, LDAF, LDB, LDX, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(IN) :: PIV(*) X INTEGER, INTENT(OUT) :: IWORK(*) X REAL(WP), INTENT(OUT) :: FERR(*), BERR(*), WORK(*) X REAL(WP), INTENT(IN) :: A(LDA,*), AF(LDAF,*), B(LDB,*) X REAL(WP), INTENT(INOUT) :: X(LDX,*) X END SUBROUTINE DGERFS X SUBROUTINE CGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, PIV, B, & X LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO ) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN) :: TRANS X INTEGER, INTENT(IN) :: LDA, LDAF, LDB, LDX, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(IN) :: PIV(*) X REAL(WP), INTENT(OUT) :: FERR(*), BERR(*), RWORK(*) X COMPLEX(WP), INTENT(OUT) :: WORK(*) X COMPLEX(WP), INTENT(IN) :: A(LDA,*), AF(LDAF,*), B(LDB,*) X COMPLEX(WP), INTENT(INOUT) :: X(LDX,*) X END SUBROUTINE CGERFS X SUBROUTINE ZGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, PIV, B, & X LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO ) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN) :: TRANS X INTEGER, INTENT(IN) :: LDA, LDAF, LDB, LDX, NRHS, N X INTEGER, INTENT(OUT) :: INFO X INTEGER, INTENT(IN) :: PIV(*) X REAL(WP), INTENT(OUT) :: FERR(*), BERR(*), RWORK(*) X COMPLEX(WP), INTENT(OUT) :: WORK(*) X COMPLEX(WP), INTENT(IN) :: A(LDA,*), AF(LDAF,*), B(LDB,*) X COMPLEX(WP), INTENT(INOUT) :: X(LDX,*) X END SUBROUTINE ZGERFS X! X SUBROUTINE SGEEQU( M,N,A,LDA,R,C,ROWCND,COLCND,AMAX,INFO ) X USE LA_PRECISION, ONLY: WP => SP X INTEGER, INTENT(IN) :: LDA, M, N X INTEGER, INTENT(OUT) :: INFO X REAL(WP), INTENT(OUT) :: AMAX, COLCND, ROWCND X REAL(WP), INTENT(IN) :: A( LDA, * ) X REAL(WP), INTENT(OUT) :: C( * ), R( * ) X END SUBROUTINE SGEEQU X SUBROUTINE DGEEQU( M,N,A,LDA,R,C,ROWCND,COLCND,AMAX,INFO ) X USE LA_PRECISION, ONLY: WP => DP X INTEGER, INTENT(IN) :: LDA, M, N X INTEGER, INTENT(OUT) :: INFO X REAL(WP), INTENT(OUT) :: AMAX, COLCND, ROWCND X REAL(WP), INTENT(IN) :: A( LDA, * ) X REAL(WP), INTENT(OUT) :: C( * ), R( * ) X END SUBROUTINE DGEEQU X SUBROUTINE CGEEQU( M,N,A,LDA,R,C,ROWCND,COLCND,AMAX,INFO ) X USE LA_PRECISION, ONLY: WP => SP X INTEGER, INTENT(IN) :: LDA, M, N X INTEGER, INTENT(OUT) :: INFO X REAL(WP), INTENT(OUT) :: AMAX, COLCND, ROWCND X COMPLEX(WP), INTENT(IN) :: A( LDA, * ) X REAL(WP), INTENT(OUT) :: C( * ), R( * ) X END SUBROUTINE CGEEQU X SUBROUTINE ZGEEQU( M,N,A,LDA,R,C,ROWCND,COLCND,AMAX,INFO ) X USE LA_PRECISION, ONLY: WP => DP X INTEGER, INTENT(IN) :: LDA, M, N X INTEGER, INTENT(OUT) :: INFO X REAL(WP), INTENT(OUT) :: AMAX, COLCND, ROWCND X COMPLEX(WP), INTENT(IN) :: A( LDA, * ) X REAL(WP), INTENT(OUT) :: C( * ), R( * ) X END SUBROUTINE ZGEEQU X! X FUNCTION SLANGE( NORM, M, N, A, LDA, WORK ) X USE LA_PRECISION, ONLY: WP => SP X REAL(WP) :: SLANGE X CHARACTER(LEN=1), INTENT(IN) :: NORM X INTEGER, INTENT(IN) :: LDA, M, N X REAL(WP), INTENT(IN) :: A( LDA, * ) X REAL(WP), INTENT(OUT) :: WORK( * ) X END FUNCTION SLANGE X FUNCTION DLANGE( NORM, M, N, A, LDA, WORK ) X USE LA_PRECISION, ONLY: WP => DP X REAL(WP) DLANGE X CHARACTER(LEN=1), INTENT(IN) :: NORM X INTEGER, INTENT(IN) :: LDA, M, N X REAL(WP), INTENT(IN) :: A( LDA, * ) X REAL(WP), INTENT(OUT) :: WORK( * ) X END FUNCTION DLANGE X FUNCTION CLANGE( NORM, M, N, A, LDA, WORK ) X USE LA_PRECISION, ONLY: WP => SP X REAL(WP) CLANGE X CHARACTER(LEN=1), INTENT(IN) :: NORM X INTEGER, INTENT(IN) :: LDA, M, N X COMPLEX(WP), INTENT(IN) :: A( LDA, * ) X REAL(WP), INTENT(OUT) :: WORK( * ) X END FUNCTION CLANGE X FUNCTION ZLANGE( NORM, M, N, A, LDA, WORK ) X USE LA_PRECISION, ONLY: WP => DP X REAL(WP) ZLANGE X CHARACTER(LEN=1), INTENT(IN) :: NORM X INTEGER, INTENT(IN) :: LDA, M, N X COMPLEX(WP), INTENT(IN) :: A( LDA, * ) X REAL(WP), INTENT(OUT) :: WORK( * ) X END FUNCTION ZLANGE X! X SUBROUTINE SGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, & X IWORK, INFO ) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN) :: NORM X INTEGER, INTENT(IN) :: LDA, N X INTEGER, INTENT(OUT) :: INFO X REAL(WP), INTENT(IN) :: ANORM X REAL(WP), INTENT(OUT) :: RCOND X INTEGER, INTENT(OUT) :: IWORK( * ) X REAL(WP), INTENT(IN) :: A( LDA, * ) X REAL(WP), INTENT(OUT) :: WORK( * ) X END SUBROUTINE SGECON X SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, & X IWORK, INFO ) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN) :: NORM X INTEGER, INTENT(IN) :: LDA, N X INTEGER, INTENT(OUT) :: INFO X REAL(WP), INTENT(IN) :: ANORM X REAL(WP), INTENT(OUT) :: RCOND X INTEGER, INTENT(OUT) :: IWORK( * ) X REAL(WP), INTENT(IN) :: A( LDA, * ) X REAL(WP), INTENT(OUT) :: WORK( * ) X END SUBROUTINE DGECON X SUBROUTINE CGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, & X RWORK, INFO ) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN) :: NORM X INTEGER, INTENT(IN) :: LDA, N X INTEGER, INTENT(OUT) :: INFO X REAL(WP), INTENT(IN) :: ANORM X REAL(WP), INTENT(OUT) :: RCOND X REAL(WP), INTENT(OUT) :: RWORK( * ) X COMPLEX(WP), INTENT(IN) :: A( LDA, * ) X COMPLEX(WP), INTENT(OUT) :: WORK( * ) X END SUBROUTINE CGECON X SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, & X RWORK, INFO ) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN) :: NORM X INTEGER, INTENT(IN) :: LDA, N X INTEGER, INTENT(OUT) :: INFO X REAL(WP), INTENT(IN) :: ANORM X REAL(WP), INTENT(OUT) :: RCOND X REAL(WP), INTENT(OUT) :: RWORK( * ) X COMPLEX(WP), INTENT(IN) :: A( LDA, * ) X COMPLEX(WP), INTENT(OUT) :: WORK( * ) X END SUBROUTINE ZGECON X! X END INTERFACE XEND MODULE LAPACK77_INTERFACES END_OF_FILE if test 16245 -ne `wc -c <'laf77mod.f90'`; then echo shar: \"'laf77mod.f90'\" unpacked with wrong size! fi # end of 'laf77mod.f90' fi if test -f 'laf90mod.f90' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'laf90mod.f90'\" else echo shar: Extracting \"'laf90mod.f90'\" \(45546 characters\) sed "s/^X//" >'laf90mod.f90' <<'END_OF_FILE' X!================================================================ X! file laf90mod.f90 X!================================================================ X! /netlib/lapack90 Preliminary version November, 1996 X! X! file: laf90mod.f90 X! for: Interface module for the LAPACK Fortran 90 routines X!================================================================ X X XMODULE LAPACK90_INTERFACES X! X INTERFACE LA_GESV X! X! Purpose X! ======= X! X! LA_GESV computes the solution to either a real or complex system of X! linear equations AX = B, X! where A is a square matrix and X and B are either rectangular X! matrices or vectors. X! X! The LU decomposition with partial pivoting and row interchanges is X! used to factor A as A = PLU, X! where P is a permutation matrix, L is unit lower triangular, and U is X! upper triangular. The factored form of A is then used to solve the X! system of equations AX = B. X! X! Arguments X! ========= X! X! SUBROUTINE LA_GESV ( A, B, IPIV, INFO ) X! (), INTENT(INOUT) :: A(:,:), X! INTEGER, INTENT(OUT), OPTIONAL :: IPIV(:) X! INTEGER, INTENT(OUT), OPTIONAL :: INFO X! where X! ::= REAL | COMPLEX X! ::= KIND(1.0) | KIND(1.0D0) X! ::= B(:,:) | B(:) X! X! ===================== X! X! A (input/output) either REAL or COMPLEX square array, X! shape (:,:), size(A,1) == size(A,2). X! On entry, the matrix A. X! On exit, the factors L and U from the factorization A = PLU; X! the unit diagonal elements of L are not stored. X! X! B (input/output) either REAL or COMPLEX rectangular array, X! shape either (:,:) or (:), size(B,1) or size(B) == size(A,1). X! On entry, the right hand side vector(s) of matrix B for the X! system of equations AX = B. X! On exit, if there is no error, the matrix of solution X! vector(s) X. X! X! IPIV Optional (output) INTEGER array, shape (:), X! size(IPIV) == size(A,1). If IPIV is present it indices that X! define the permutation matrix P; row i of the matrix was X! interchanged with row IPIV(i). X! X! INFO Optional (output) INTEGER. X! If INFO is present X! = 0: successful exit X! < 0: if INFO = -k, the k-th argument had an illegal value X! > 0: if INFO = k, U(k,k) is exactly zero. The factorization X! has been completed, but the factor U is exactly X! singular, so the solution could not be computed. X! If INFO is not present and an error occurs, then the program is X! terminated with an error message. X! X! ===================================================================== X X SUBROUTINE SGESV_F90(A,B,IPIV,INFO) X USE LA_PRECISION, ONLY: SP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(OUT), OPTIONAL :: IPIV(:) X REAL(SP), INTENT(INOUT) :: A(:,:), B(:,:) X END SUBROUTINE SGESV_F90 X SUBROUTINE S1GESV_F90(A,B,IPIV,INFO) X USE LA_PRECISION, ONLY: SP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(OUT), OPTIONAL :: IPIV(:) X REAL(SP), INTENT(INOUT) :: A(:,:), B(:) X END SUBROUTINE S1GESV_F90 X SUBROUTINE DGESV_F90(A,B,IPIV,INFO) X USE LA_PRECISION, ONLY: DP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(OUT), OPTIONAL :: IPIV(:) X REAL(DP), INTENT(INOUT) :: A(:,:), B(:,:) X END SUBROUTINE DGESV_F90 X SUBROUTINE D1GESV_F90(A,B,IPIV,INFO) X USE LA_PRECISION, ONLY: DP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(OUT), OPTIONAL :: IPIV(:) X REAL(DP), INTENT(INOUT) :: A(:,:), B(:) X END SUBROUTINE D1GESV_F90 X SUBROUTINE CGESV_F90(A,B,IPIV,INFO) X USE LA_PRECISION, ONLY: SP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(OUT), OPTIONAL :: IPIV(:) X COMPLEX(SP), INTENT(INOUT) :: A(:,:), B(:,:) X END SUBROUTINE CGESV_F90 X SUBROUTINE C1GESV_F90(A,B,IPIV,INFO) X USE LA_PRECISION, ONLY: SP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(OUT), OPTIONAL :: IPIV(:) X COMPLEX(SP), INTENT(INOUT) :: A(:,:), B(:) X END SUBROUTINE C1GESV_F90 X SUBROUTINE ZGESV_F90(A,B,IPIV,INFO) X USE LA_PRECISION, ONLY: DP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(OUT), OPTIONAL :: IPIV(:) X COMPLEX(DP), INTENT(INOUT) :: A(:,:), B(:,:) X END SUBROUTINE ZGESV_F90 X SUBROUTINE Z1GESV_F90(A,B,IPIV,INFO) X USE LA_PRECISION, ONLY: DP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(OUT), OPTIONAL :: IPIV(:) X COMPLEX(DP), INTENT(INOUT) :: A(:,:), B(:) X END SUBROUTINE Z1GESV_F90 X END INTERFACE X! X INTERFACE LA_GESVX X! X! Purpose X! ======= X! X! LA_GESVX computes the solution to a either real or complex system of X! linear equations AX = B, X! where A is a square matrix and X and B are either rectangular X! matrices or vectors. X! LA_GESVX is an expert driver routine, which can also perform the X! following functions: X! - solve A^T X = B or A^H X = B, X! - estimate the condition number of A, X! - return the pivot growth factor, X! - refine the solution and computes forward and backward error X! bounds, X! - equilibrate the system if A is poorly scaled. X! X! Arguments X! ========= X! SUBROUTINE LA_GESVX (A, B, X, AF, IPIV, FACT, TRANS, EQUED, R, C, X! FERR, BERR, RCOND, RPVGRW, INFO) X! (), INTENT(INOUT) :: A(:,:), X! (), INTENT(OUT) :: X! (), INTENT(INOUT), OPTIONAL :: AF(:,:) X! INTEGER, INTENT(INOUT), OPTIONAL :: IPIV(:) X! CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED X! REAL(), INTENT(INOUT), OPTIONAL :: R(:), C(:) X! REAL(), INTENT(OUT), OPTIONAL :: , RCOND, RPVGRW X! CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: FACT, TRANS X! INTEGER, INTENT(OUT), OPTIONAL :: INFO X! where X! ::= REAL | COMPLEX X! ::= KIND(1.0) | KIND(1.0D0) X! ::= B(:,:) | B(:) X! ::= X(:,:) | X(:) X! ::= FERR(:), BERR(:) | FERR, BERR X! X! Description X! =========== X! X! The following steps are performed: X! X! 1. If FACT is not present or FACT = 'N', and EQUED is present, X! real scaling factors are computed to equilibrate the system: X! TRANS = 'N': diag(R) A diag(C) inv(diag(C)) X = diag(R) B X! TRANS = 'T': (diag(R) A diag(C))^T inv(diag(R)) X = diag(C) B X! TRANS = 'C': (diag(R) A diag(C))^H inv(diag(R)) X = diag(C) B X! Whether or not the system will be equilibrated depends on the X! scaling of the matrix A, but if equilibration is used, A is X! overwritten by diag(R) A diag(C) and B by diag(R) B (if TRANS='N') X! or diag(C) B (if TRANS = 'T' or 'C'). X! X! 2. If FACT = 'N', the LU decomposition is used to factor the X! matrix A (after equilibration if EQUED is present) as A = PLU, X! where P is a permutation matrix, L is a unit lower triangular X! matrix, and U is upper triangular. X! X! 3. The factored form of A is used to estimate the condition number X! of the matrix A. If the reciprocal of the condition number is X! less than machine precision, steps 4-6 are skipped. X! X! 4. The system of equations is solved for X using the factored form X! of A. X! X! 5. Iterative refinement is applied to improve the computed solution X! matrix and calculate error bounds and backward error estimates X! for it. X! X! 6. If equilibration was used, the matrix X is premultiplied by X! diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so X! that it solves the original system before equilibration. X! X! Arguments X! ========= X! X! A (input/output) either REAL or COMPLEX square array, X! shape (:,:), size(A,1) == size(A,2). X! If FACT is not present or FACT = 'N', X! On entry, the matrix A. X! On exit, if EQUED is present, the matrix A may have been X! overwritten by the equilibrated matrix (see EQUED). X! If FACT is present and FACT = 'F', X! On entry, the matrix A, possibly equilibrated in a previous X! call to LA_GESVX (see EQUED). X! On exit, A is unchanged. X! X! B (input/output) either REAL or COMPLEX rectangular array, X! shape either (:,:) or (:), size(B,1) or size(B) == size(A,1). X! On entry, the right hand side vector(s) of matrix B for the X! system of equations AX = B. X! On exit, if EQUED is present, B may have been scaled in X! accordance with the equilibration of A (see EQUED); X! otherwise, B is unchanged. X! X! X (output) either REAL or COMPLEX rectangular array, X! shape either (:,:) or (:), size(X,1) or size(X) == size(A,1). X! If INFO = 0, the solution matrix X to the original X! system of equations. Note that X always returns the solution X! to the original system of equations; if equilibration has been X! performed (EQUED is present and EQUED /= 'N'), this does not X! correspond to the scaled A and B. X! X! AF Optional (input/output) either REAL or COMPLEX square array, X! shape (:,:), size(AF,1) == size(AF,2) == size(A,1). X! If FACT is not present or FACT = 'N', then AF is an output X! argument and returns the factors L and U from the factorization X! A = PLU of the original matrix A, possibly equilibrated if X! EQUED is present. X! If FACT is present and FACT = 'F', then AF is an input argument X! (and must be present); on entry, it must contain the factors X! L and U of A (possibly equilibrated if EQUED is present), X! returned by a previous call to LA_GESVX. X! X! IPIV Optional (input/output) INTEGER array, shape (:), X! size(IPIV) == size(A,1). X! If FACT is not present or FACT = 'N', then IPIV is an output X! argument and returns the pivot indices from the factorization X! A = PLU of the original matrix A, possibly equilibrated if X! EQUED is present. X! If FACT is present and FACT = 'F', then IPIV is an input argument X! (and must be present); on entry, it must contain the pivot X! indices from the factorization of A (possibly equilibrated X! if EQUED is present), returned by a previous call to LA_GESVX. X! X! TRANS Optional (input) CHARACTER*1 X! If TRANS is present, it specifies the form of the system X! of equations: X! = 'N': A X = B (No transpose) X! = 'T': A^T X = B (Transpose) X! = 'C': A^H X = B (Conjugate transpose = Transpose) X! otherwise TRANS = 'N' is assumed. X! X! FACT Optional (input) CHARACTER*1 X! Specifies whether or not the factored form of the matrix A is X! supplied on entry. X! If FACT is present then: X! = 'N': the matrix A will be equilibrated if EQUED is present, X! then copied to AF and factored. X! = 'F': on entry, AF and IPIV must contain the factored form X! of A (possibly equilibrated if EQUED is present). X! otherwise FACT = 'N' is assumed. X! X! EQUED Optional (input/output) CHARACTER*1 X! If FACT is not present or FACT = 'N', then EQUED is an output X! argument. If it is present, then the matrix is equilibrated, X! and on exit EQUED specifies the scaling of A which has X! actually been performed: X! = 'N': No equilibration. X! = 'R': Row equilibration, i.e., A has been premultiplied X! by diag(R); also B has been premultiplied by diag(R) X! if TRANS = 'N'. X! = 'C': Column equilibration, i.e., A has been postmultiplied X! by diag(C); also B has been postmultiplied by diag(C) X! if TRANS = 'T' or 'C'. X! = 'B': Both row and column equilibration: combines the X! effects of EQUED = 'R' and EQUED = 'C'. X! If FACT is present and FACT = 'F', then EQUED is an input X! argument; if it is present, it specifies the equilibration X! of A which was performed in a previous call to LA_GESVX with X! FACT not present or FACT = 'N'. X! X! R Optional (input/output) REAL array, shape (:), X! size(R) == size(A,1). X! R must be present if EQUED is present and EQUED = 'R' or 'B'; X! R is not referenced if EQUED = 'N' or 'C'. X! If FACT is not present or FACT = 'N', then R is an output X! argument. If EQUED = 'R' or 'B', R returns the row scale X! factors for equilibrating A. X! If FACT is present and FACT = 'F', then R is an input X! argument. If EQUED = 'R' or 'B', R must contain the row scale X! factors for equilibrating A, returned by a previous call to X! LA_GESVX; each element of R must be positive. X! X! C Optional (input/output) REAL array, shape (:), X! size(C) == size(A,1). X! C must be present if EQUED is present and EQUED = 'C' or 'B'; X! C is not referenced if EQUED = 'N' or 'R'. X! If FACT is not present or FACT = 'N', then C is an output X! argument. If EQUED = 'C' or 'B', C returns the column scale X! factors for equilibrating A. X! If FACT is present and FACT = 'F', then C is an input X! argument. If EQUED = 'C' or 'B', C must contain the column X! scale factors for equilibrating A, returned by a previous X! call to LA_GESVX; each element of C must be positive. X! X! FERR Optional (output) either REAL array of shape (:) or REAL X! scalar. If it is an array, size(FERR) == size(X,2). X! The estimated forward error bound for each solution vector X! X(j) (the j-th column of the solution matrix X). X! If XTRUE is the true solution corresponding to X(j), FERR(j) X! is an estimated upper bound for the magnitude of the largest X! element in (X(j) - XTRUE) divided by the magnitude of the X! largest element in X(j). The estimate is as reliable as X! the estimate for RCOND, and is almost always a slight X! overestimate of the true error. X! X! BERR Optional (output) either REAL array of shape (:) or REAL X! scalar. If it is an array, size(BERR) == size(X,2). X! The componentwise relative backward error of each solution X! vector X(j) (i.e., the smallest relative change in X! any element of A or B that makes X(j) an exact solution). X! X! RCOND Optional (output) REAL X! The estimate of the reciprocal condition number of the matrix X! A after equilibration (if done). If RCOND is less than the X! machine precision (in particular, if RCOND = 0), the matrix X! is singular to working precision. This condition is X! indicated by a return code of INFO > 0, and the solution and X! error bounds are not computed. X! X! RPVGRW Optional (output) REAL. X! The reciprocal pivot growth factor norm(A)/norm(U). X! If RPVGRW is much less than 1, then the stability X! of the LU factorization of the (equilibrated) matrix A X! could be poor. This also means that the solution X, condition X! estimator RCOND, and forward error bound FERR could be X! unreliable. If factorization fails with 0 0: if INFO = k, and k is X! <= N: U(k,k) is exactly zero. The factorization has X! been completed, but the factor U is exactly X! singular, so the solution and error bounds X! could not be computed. X! = N+1: RCOND is less than machine precision. The X! factorization has been completed, but the X! matrix is singular to working precision, and X! the solution and error bounds have not been X! computed. X! If INFO is not present and an error occurs, then the program is X! terminated with an error message. X! X! ===================================================================== X! X SUBROUTINE SGESVX_F90(A, B, X, AF, IPIV, FACT, TRANS, & X EQUED, R, C, FERR, BERR, RCOND, RPVGRW, INFO) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS, FACT X CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(INOUT), OPTIONAL :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: RCOND, RPVGRW X REAL(WP), INTENT(OUT), OPTIONAL :: FERR(:), BERR(:) X REAL(WP), INTENT(INOUT), OPTIONAL :: AF(:,:), C(:), R(:) X REAL(WP), INTENT(INOUT) :: A(:,:), B(:,:) X REAL(WP), INTENT(OUT) :: X(:,:) X END SUBROUTINE SGESVX_F90 X SUBROUTINE S1GESVX_F90(A, B, X, AF, IPIV, FACT, TRANS, & X EQUED, R, C, FERR, BERR, RCOND, RPVGRW, INFO) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS, FACT X CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(INOUT), OPTIONAL :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: RCOND, FERR, BERR, RPVGRW X REAL(WP), INTENT(INOUT), OPTIONAL :: AF(:,:), C(:), R(:) X REAL(WP), INTENT(INOUT) :: A(:,:), B(:) X REAL(WP), INTENT(OUT) :: X(:) X END SUBROUTINE S1GESVX_F90 X SUBROUTINE DGESVX_F90(A, B, X, AF, IPIV, FACT, TRANS, & X EQUED, R, C, FERR, BERR, RCOND, RPVGRW, INFO) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS, FACT X CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(INOUT), OPTIONAL :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: RCOND, RPVGRW X REAL(WP), INTENT(OUT), OPTIONAL :: FERR(:), BERR(:) X REAL(WP), INTENT(INOUT), OPTIONAL :: AF(:,:), C(:), R(:) X REAL(WP), INTENT(INOUT) :: A(:,:), B(:,:) X REAL(WP), INTENT(OUT) :: X(:,:) X END SUBROUTINE DGESVX_F90 X SUBROUTINE D1GESVX_F90(A, B, X, AF, IPIV, FACT, TRANS, & X EQUED, R, C, FERR, BERR, RCOND, RPVGRW, INFO) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS, FACT X CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(INOUT), OPTIONAL :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: RCOND, FERR, BERR, RPVGRW X REAL(WP), INTENT(INOUT), OPTIONAL :: AF(:,:), C(:), R(:) X REAL(WP), INTENT(INOUT) :: A(:,:), B(:) X REAL(WP), INTENT(OUT) :: X(:) X END SUBROUTINE D1GESVX_F90 X SUBROUTINE CGESVX_F90(A, B, X, AF, IPIV, FACT, TRANS, & X EQUED, R, C, FERR, BERR, RCOND, RPVGRW, INFO) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS, FACT X CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(INOUT), OPTIONAL :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: RCOND, RPVGRW X REAL(WP), INTENT(OUT), OPTIONAL :: FERR(:), BERR(:) X REAL(WP), INTENT(INOUT), OPTIONAL :: C(:), R(:) X COMPLEX(WP), INTENT(INOUT), OPTIONAL :: AF(:,:) X COMPLEX(WP), INTENT(INOUT) :: A(:,:), B(:,:) X COMPLEX(WP), INTENT(OUT) :: X(:,:) X END SUBROUTINE CGESVX_F90 X SUBROUTINE C1GESVX_F90(A, B, X, AF, IPIV, FACT, TRANS, & X EQUED, R, C, FERR, BERR, RCOND, RPVGRW, INFO) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS, FACT X CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(INOUT), OPTIONAL :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: RCOND, FERR, BERR, RPVGRW X REAL(WP), INTENT(INOUT), OPTIONAL :: C(:), R(:) X COMPLEX(WP), INTENT(INOUT), OPTIONAL :: AF(:,:) X COMPLEX(WP), INTENT(INOUT) :: A(:,:), B(:) X COMPLEX(WP), INTENT(OUT) :: X(:) X END SUBROUTINE C1GESVX_F90 X SUBROUTINE ZGESVX_F90(A, B, X, AF, IPIV, FACT, TRANS, & X EQUED, R, C, FERR, BERR, RCOND, RPVGRW, INFO) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS, FACT X CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(INOUT), OPTIONAL :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: RCOND, RPVGRW X REAL(WP), INTENT(OUT), OPTIONAL :: FERR(:), BERR(:) X REAL(WP), INTENT(INOUT), OPTIONAL :: C(:), R(:) X COMPLEX(WP), INTENT(INOUT), OPTIONAL :: AF(:,:) X COMPLEX(WP), INTENT(INOUT) :: A(:,:), B(:,:) X COMPLEX(WP), INTENT(OUT) :: X(:,:) X END SUBROUTINE ZGESVX_F90 X SUBROUTINE Z1GESVX_F90(A, B, X, AF, IPIV, FACT, TRANS, & X EQUED, R, C, FERR, BERR, RCOND, RPVGRW, INFO) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS, FACT X CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(INOUT), OPTIONAL :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: RCOND, FERR, BERR, RPVGRW X REAL(WP), INTENT(INOUT), OPTIONAL :: C(:), R(:) X COMPLEX(WP), INTENT(INOUT), OPTIONAL :: AF(:,:) X COMPLEX(WP), INTENT(INOUT) :: A(:,:), B(:) X COMPLEX(WP), INTENT(OUT) :: X(:) X END SUBROUTINE Z1GESVX_F90 X END INTERFACE X! X INTERFACE LA_GETRF X! X! Purpose X! ======= X! X! LA_GETRF computes an LU factorization of a general regtangle X! matrix A using partial pivoting with row interchanges. X! X! The factorization has the form A = PLU X! where P is a permutation matrix, L is lower triangular X! with unit diagonal elements (lower trapezoidal if m > n), X! and U is upper triangular (upper trapezoidal if m < n), X! where m=size(A,1) and n=size(A,2). X! X! When A is square (m = n), LA_GETRF optionally estimates the X! reciprocal of the condition number of a general matrix A, in either X! the 1-norm or the infinity-norm. X! An estimate is obtained for norm(inv(A)), and the reciprocal of the X! condition number is computed as RCOND = 1 / (norm(A) * norm(inv(A))). X! X! Arguments X! ========= X! X! SUBROUTINE LA_GETRF ( A, IPIV, RCOND, NORM, INFO ) X! (), INTENT(INOUT) :: A(:,:) X! INTEGER, INTENT( OUT ) :: IPIV( : ) X! REAL(), INTENT( OUT ), OPTIONAL :: RCOND X! CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: NORM X! INTEGER, INTENT(OUT), OPTIONAL :: INFO X! where X! ::= REAL | COMPLEX X! ::= KIND(1.0) | KIND(1.0D0) X! X! ==================== X! X! A (input/output) either REAL or COMPLEX array, shape (:,:). X! On entry, the matrix A. X! On exit, the factors L and U from the factorization A = PLU; X! the unit diagonal elements of L are not stored. X! X! IPIV (output) INTEGER array, shape (:), X! size(IPIV) == min(size(A,1),size(A,2)). X! IPIV indices that define the permutation matrix P; X! row i of the matrix A was interchanged with row IPIV(i). X! X! RCOND Optional (output) REAL X! The reciprocal of the condition number of the matrix A for X! the case m = n, computed as RCOND = 1/(norm(A) * norm(inv(A))). X! RCOND should be present if NORM is present. X! If m /= n then RCOND is returned as zero. X! X! NORM Optional (input) CHARACTER*1 X! Specifies whether the 1-norm condition number or the X! infinity-norm condition number is required: X! = '1', 'O' or 'o': 1-norm; X! = 'I', 'i': infinity-norm. X! If NORM is not present, the 1-norm is used. X! X! INFO Optional (output) INTEGER. X! If INFO is present X! = 0: successful exit X! < 0: if INFO = -k, the k-th argument had an illegal value X! > 0: if INFO = k, U(k,k) is exactly zero. The factorization X! has been completed, but the factor U is exactly X! singular, so the solution could not be computed. X! If INFO is not present and an error occurs, then the program is X! terminated with an error message. X! X! ===================================================================== X! X SUBROUTINE SGETRF_F90( A,IPIV,RCOND,NORM,INFO ) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: NORM X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT( OUT ) :: IPIV( : ) X REAL(WP), INTENT( INOUT ) :: A( :, : ) X REAL(WP), INTENT( OUT ), OPTIONAL :: RCOND X END SUBROUTINE SGETRF_F90 X SUBROUTINE DGETRF_F90( A,IPIV,RCOND,NORM,INFO ) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: NORM X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT( OUT ) :: IPIV( : ) X REAL(WP), INTENT( INOUT ) :: A( :, : ) X REAL(WP), INTENT( OUT ), OPTIONAL :: RCOND X END SUBROUTINE DGETRF_F90 X SUBROUTINE CGETRF_F90( A,IPIV,RCOND,NORM,INFO ) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: NORM X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT( OUT ) :: IPIV( : ) X COMPLEX(WP), INTENT( INOUT ) :: A( :, : ) X REAL(WP), INTENT( OUT ), OPTIONAL :: RCOND X END SUBROUTINE CGETRF_F90 X SUBROUTINE ZGETRF_F90( A,IPIV,RCOND,NORM,INFO ) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: NORM X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT( OUT ) :: IPIV( : ) X COMPLEX(WP), INTENT( INOUT ) :: A( :, : ) X REAL(WP), INTENT( OUT ), OPTIONAL :: RCOND X END SUBROUTINE ZGETRF_F90 X END INTERFACE X! X INTERFACE LA_GETRS X! X! Purpose X! ======= X! X! LA_GETRS solves a system of linear equations X! A X = B, A^T X = B or A^H X = B X! with a general square matrix A using the LU factorization computed X! by LA_GETRF. X! X! Arguments X! ========= X! SUBROUTINE LA_GETRS (A, IPIV, B, TRANS, INFO) X! (), INTENT(IN) :: A(:,:) X! (), INTENT(INOUT) :: X! INTEGER, INTENT(IN) :: IPIV(:) X! CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X! INTEGER, INTENT(OUT), OPTIONAL :: INFO X! where X! ::= REAL | COMPLEX X! ::= KIND(1.0) | KIND(1.0D0) X! ::= B(:,:) | B(:) X! X! ===================== X! X! A (input) either REAL or COMPLEX square array, X! shape (:,:), size(A,1) == size(A,2). X! The factors L and U from the factorization A = PLU as computed X! by LA_GETRF. X! X! IPIV (input) INTEGER array, shape (:), size(IPIV) == size(A,1). X! The pivot indices from LA_GETRF; for 1<=i<=size(A,1), row i X! of the matrix was interchanged with row IPIV(i). X! X! B (input/output) either REAL or COMPLEX rectangular array, X! shape either (:,:) or (:), size(B,1) or size(B) == size(A,1). X! On entry, the right hand side vector(s) of matrix B for X! the system of equations AX = B. X! On exit, if there is no error, the matrix of solution X! vector(s) X. X! X! TRANS Optional (input) CHARACTER*1 X! If TRANS is present, it specifies the form of the system X! of equations: X! = 'N': A X = B (No transpose) X! = 'T': A^T X = B (Transpose) X! = 'C': A^H X = B (Conjugate transpose = Transpose) X! otherwise TRANS = 'N' is assumed. X! X! INFO Optional (output) INTEGER. X! If INFO is present X! = 0: successful exit X! < 0: if INFO = -k, the k-th argument had an illegal value X! If INFO is not present and an error occurs, then the program is X! terminated with an error message. X! X! ===================================================================== X! X SUBROUTINE SGETRS_F90(A,IPIV,B,TRANS,INFO) X USE LA_PRECISION, ONLY: SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(SP), INTENT(IN) :: A(:,:) X REAL(SP), INTENT(INOUT) :: B(:,:) X END SUBROUTINE SGETRS_F90 X SUBROUTINE DGETRS_F90(A,IPIV,B,TRANS,INFO) X USE LA_PRECISION, ONLY: DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(DP), INTENT(IN) :: A(:,:) X REAL(DP), INTENT(INOUT) :: B(:,:) X END SUBROUTINE DGETRS_F90 X SUBROUTINE CGETRS_F90(A,IPIV,B,TRANS,INFO) X USE LA_PRECISION, ONLY: SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X COMPLEX(SP), INTENT(IN) :: A(:,:) X COMPLEX(SP), INTENT(INOUT) :: B(:,:) X END SUBROUTINE CGETRS_F90 X SUBROUTINE ZGETRS_F90(A,IPIV,B,TRANS,INFO) X USE LA_PRECISION, ONLY: DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X COMPLEX(DP), INTENT(IN) :: A(:,:) X COMPLEX(DP), INTENT(INOUT) :: B(:,:) X END SUBROUTINE ZGETRS_F90 X SUBROUTINE S1GETRS_F90(A,IPIV,B,TRANS,INFO) X USE LA_PRECISION, ONLY: SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(SP), INTENT(IN) :: A(:,:) X REAL(SP), INTENT(INOUT) :: B(:) X END SUBROUTINE S1GETRS_F90 X SUBROUTINE D1GETRS_F90(A,IPIV,B,TRANS,INFO) X USE LA_PRECISION, ONLY: DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(DP), INTENT(IN) :: A(:,:) X REAL(DP), INTENT(INOUT) :: B(:) X END SUBROUTINE D1GETRS_F90 X SUBROUTINE C1GETRS_F90(A,IPIV,B,TRANS,INFO) X USE LA_PRECISION, ONLY: SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X COMPLEX(SP), INTENT(IN) :: A(:,:) X COMPLEX(SP), INTENT(INOUT) :: B(:) X END SUBROUTINE C1GETRS_F90 X SUBROUTINE Z1GETRS_F90(A,IPIV,B,TRANS,INFO) X USE LA_PRECISION, ONLY: DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X COMPLEX(DP), INTENT(IN) :: A(:,:) X COMPLEX(DP), INTENT(INOUT) :: B(:) X END SUBROUTINE Z1GETRS_F90 X END INTERFACE X! X INTERFACE LA_GETRI X! X! Purpose X! ======= X! X! LA_GETRI computes the inverse of a matrix using the LU factorization X! computed by LA_GETRF. X! X! Arguments X! ========= X! SUBROUTINE LA_GETRI (A, IPIV, INFO) X! (), INTENT(INOUT) :: A(:,:) X! INTEGER, INTENT(IN) :: IPIV(:) X! INTEGER, INTENT(OUT), OPTIONAL :: INFO X! where X! ::= REAL | COMPLEX X! ::= KIND(1.0) | KIND(1.0D0) X! X! ===================== X! X! A (input/output) either REAL or COMPLEX square array, shape (:,:), X! size(A,1) == size(A,2). X! On entry contains the factors L and U from the factorization X! A = PLU as computed by LA_GETRF. X! On exit, if INFO = 0, the inverse of the original matrix A. X! X! IPIV (input) INTEGER array, shape (:), size(IPIV) == size(A,1). X! The pivot indices from LA_GETRF; for 1<=i<=size(A,1), row i of X! the matrix was interchanged with row IPIV(i). X! X! INFO Optional (output) INTEGER. X! If INFO is present X! = 0: successful exit X! < 0: if INFO = -k, the k-th argument had an illegal value X! > 0: if INFO = k, U(k,k) is exactly zero. The matrix is X! singular and its inverse could not be computed. X! If INFO is not present and an error occurs, then the program is X! terminated with an error message. X! X! ===================================================================== X! X SUBROUTINE SGETRI_F90(A,IPIV,INFO) X USE LA_PRECISION, ONLY: WP => SP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(WP), INTENT(INOUT) :: A(:,:) X END SUBROUTINE SGETRI_F90 X SUBROUTINE DGETRI_F90(A,IPIV,INFO) X USE LA_PRECISION, ONLY: WP => DP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(WP), INTENT(INOUT) :: A(:,:) X END SUBROUTINE DGETRI_F90 X SUBROUTINE CGETRI_F90(A,IPIV,INFO) X USE LA_PRECISION, ONLY: WP => SP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X COMPLEX(WP), INTENT(INOUT) :: A(:,:) X END SUBROUTINE CGETRI_F90 X SUBROUTINE ZGETRI_F90(A,IPIV,INFO) X USE LA_PRECISION, ONLY: WP => DP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X COMPLEX(WP), INTENT(INOUT) :: A(:,:) X END SUBROUTINE ZGETRI_F90 X END INTERFACE X! X INTERFACE LA_GERFS X! X! Purpose X! ======= X! X! LA_GERFS improves the computed solution X of a system of linear X! equations A X = B or A^T X = B X! and provides error bounds and backward error estimates for X! the solution. LA_GERFS uses the LU factors computed by LA_GETRF. X! X! Arguments X! ========= X! SUBROUTINE LA_GERFS (A, AF, IPIV, B, X, TRANS, FERR, BERR, INFO) X! (), INTENT(IN) :: A(:,:), AF(:,:), X! INTEGER, INTENT(IN) :: IPIV(:) X! (), INTENT(INOUT) :: X! REAL(), INTENT(OUT), OPTIONAL :: X! CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X! INTEGER, INTENT(OUT), OPTIONAL :: INFO X! where X! ::= REAL | COMPLEX X! ::= KIND(1.0) | KIND(1.0D0) X! ::= B(:,:) | B(:) X! ::= X(:,:) | X(:) X! ::= FERR(:), BERR(:) | FERR, BERR X! X! ===================== X! X! A (input) either REAL or COMPLEX square array, X! shape (:,:), size(A,1) == size(A,2). X! The original matrix A. X! X! AF (input) either REAL or COMPLEX square array, X! shape (:,:), size(AF,1) == size(AF,2) == size(A,1). X! The factors L and U from the factorization A = PLU X! as computed by LA_GETRF. X! X! IPIV (input) INTEGER array, shape (:), size(IPIV) == size(A,1). X! The pivot indices from LA_GETRF; for 1<=i<=size(A,1), row i X! of the matrix was interchanged with row IPIV(i). X! X! B (input) either REAL or COMPLEX rectangular array, X! shape either (:,:) or (:), size(B,1) or size(B) == size(A,1). X! The right hand side vector(s) of matrix B for X! the system of equations AX = B. X! X! X (input/output) either REAL or COMPLEX rectangular array, X! shape either (:,:) or (:), size(X,1) or size(X) == size(A,1). X! On entry, the solution matrix X, as computed by LA_GETRS. X! On exit, the improved solution matrix X. X! X! TRANS Optional (input) CHARACTER*1 X! If TRANS is present, it specifies the form of the system X! of equations: X! = 'N': A X = B (No transpose) X! = 'T': A^T X = B (Transpose) X! = 'C': A^H X = B (Conjugate transpose = Transpose) X! otherwise TRANS = 'N' is assumed. X! X! FERR Optional (output) either REAL array of shape (:) or REAL X! scalar. If it is an array, size(FERR) == size(X,2). X! The estimated forward error bound for each solution vector X! X(j) (the j-th column of the solution matrix X). X! If XTRUE is the true solution corresponding to X(j), FERR(j) X! is an estimated upper bound for the magnitude of the largest X! element in (X(j) - XTRUE) divided by the magnitude of the X! largest element in X(j). The estimate is as reliable as X! the estimate for RCOND, and is almost always a slight X! overestimate of the true error. X! X! BERR Optional (output) either REAL array of shape (:) or REAL X! scalar. If it is an array, size(BERR) == size(X,2). X! The componentwise relative backward error of each solution X! vector X(j) (i.e., the smallest relative change in X! any element of A or B that makes X(j) an exact solution). X! X! INFO Optional (output) INTEGER. X! If INFO is present X! = 0: successful exit X! < 0: if INFO = -k, the k-th argument had an illegal value X! If INFO is not present and an error occurs, then the program is X! terminated with an error message. X! X! Internal Parameters X! =================== X! X! ITMAX is the maximum number of steps of iterative refinement. X! It is set to 5 in the LAPACK77 subroutines X! X! ===================================================================== X! X SUBROUTINE SGERFS_F90(A,AF,IPIV,B,X,TRANS,FERR,BERR,INFO) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: FERR(:), BERR(:) X REAL(WP), INTENT(IN) :: A(:,:), AF(:,:), B(:,:) X REAL(WP), INTENT(INOUT) :: X(:,:) X END SUBROUTINE SGERFS_F90 X SUBROUTINE S1GERFS_F90(A,AF,IPIV,B,X,TRANS,FERR,BERR,INFO) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: FERR, BERR X REAL(WP), INTENT(IN) :: A(:,:), AF(:,:), B(:) X REAL(WP), INTENT(INOUT) :: X(:) X END SUBROUTINE S1GERFS_F90 X SUBROUTINE DGERFS_F90(A,AF,IPIV,B,X,TRANS,FERR,BERR,INFO) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: FERR(:), BERR(:) X REAL(WP), INTENT(IN) :: A(:,:), AF(:,:), B(:,:) X REAL(WP), INTENT(INOUT) :: X(:,:) X END SUBROUTINE DGERFS_F90 X SUBROUTINE D1GERFS_F90(A,AF,IPIV,B,X,TRANS,FERR,BERR,INFO) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: FERR, BERR X REAL(WP), INTENT(IN) :: A(:,:), AF(:,:), B(:) X REAL(WP), INTENT(INOUT) :: X(:) X END SUBROUTINE D1GERFS_F90 X SUBROUTINE CGERFS_F90(A,AF,IPIV,B,X,TRANS,FERR,BERR,INFO) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: FERR(:), BERR(:) X COMPLEX(WP), INTENT(IN) :: A(:,:), AF(:,:), B(:,:) X COMPLEX(WP), INTENT(INOUT) :: X(:,:) X END SUBROUTINE CGERFS_F90 X SUBROUTINE C1GERFS_F90(A,AF,IPIV,B,X,TRANS,FERR,BERR,INFO) X USE LA_PRECISION, ONLY: WP => SP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: FERR, BERR X COMPLEX(WP), INTENT(IN) :: A(:,:), AF(:,:), B(:) X COMPLEX(WP), INTENT(INOUT) :: X(:) X END SUBROUTINE C1GERFS_F90 X SUBROUTINE ZGERFS_F90(A,AF,IPIV,B,X,TRANS,FERR,BERR,INFO) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: FERR(:), BERR(:) X COMPLEX(WP), INTENT(IN) :: A(:,:), AF(:,:), B(:,:) X COMPLEX(WP), INTENT(INOUT) :: X(:,:) X END SUBROUTINE ZGERFS_F90 X SUBROUTINE Z1GERFS_F90(A,AF,IPIV,B,X,TRANS,FERR,BERR,INFO) X USE LA_PRECISION, ONLY: WP => DP X CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: TRANS X INTEGER, INTENT(OUT), OPTIONAL :: INFO X INTEGER, INTENT(IN) :: IPIV(:) X REAL(WP), INTENT(OUT), OPTIONAL :: FERR, BERR X COMPLEX(WP), INTENT(IN) :: A(:,:), AF(:,:), B(:) X COMPLEX(WP), INTENT(INOUT) :: X(:) X END SUBROUTINE Z1GERFS_F90 X END INTERFACE X! X INTERFACE LA_GEEQU X! X! Purpose X! ======= X! X! LA_GEEQU computes row and column scalings intended to equilibrate a X! rectangle matrix A and reduce its condition number. R returns the X! row scale factors and C the column scale factors, chosen to try to X! make the largest entry in each row and column of the matrix B with X! elements B(i,j) = R(i) A(i,j) C(j) have absolute value 1. X! X! R(i) and C(j) are restricted to be between SMLNUM = smallest safe X! number and BIGNUM = largest safe number. Use of these scaling X! factors is not guaranteed to reduce the condition number of A but X! works well in practice. X! X! Arguments X! ========= X! X! SUBROUTINE LA_GEEQU ( A, R, C, ROWCND, COLCND, AMAX, INFO ) X! (), INTENT(IN) :: A(:,:) X! REAL(), INTENT( OUT ) :: R(:), C(:) X! REAL(), INTENT( OUT ), OPTIONAL :: ROWCND, COLCND, AMAX X! INTEGER, INTENT(OUT), OPTIONAL :: INFO X! where X! ::= REAL | COMPLEX X! ::= KIND(1.0) | KIND(1.0D0) X! X! ==================== X! X! A (input) either REAL or COMPLEX array, shape (:,:). X! The matrix A, whose equilibration factors are to be computed. X! X! R (output) REAL array, shape (:), size(R) == size(A,1). X! If INFO = 0 or INFO > size(A,1), R contains the row X! scale factors for A. X! X! C (output) REAL array, shape (:), size(C) == size(A,2). X! If INFO = 0, C contains the column scale factors for A. X! X! ROWCND Optional (output) REAL. X! If INFO = 0 or INFO > size(A,1), ROWCND contains the ratio X! of the smallest R(i) to the largest R(i). If ROWCND >= 0.1 X! and AMAX is neither too large nor too small, it is not worth X! scaling by R. X! X! COLCND Optional (output) REAL. X! If INFO = 0, COLCND contains the ratio of the smallest X! C(i) to the largest C(i). If COLCND >= 0.1, it is not X! worth scaling by C. X! X! AMAX Optional (output) REAL. X! Absolute value of largest matrix element. If AMAX is very X! close to overflow or very close to underflow, the matrix X! should be scaled. X! X! INFO Optional (output) INTEGER X! If INFO is present X! = 0: successful exit X! < 0: if INFO = -k, the k-th argument had an illegal value X! > 0: if INFO = k, and k is X! <= M: the k-th row of A is exactly zero X! > M: the (k-M)-th column of A is exactly zero X! where M = size(A,1) X! If INFO is not present and an error occurs, then the program is X! terminated with an error message. X! X! ===================================================================== X SUBROUTINE SGEEQU_F90( A, R, C, ROWCND, COLCND, AMAX, INFO ) X USE LA_PRECISION, ONLY: WP => SP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X REAL(WP), INTENT( IN ) :: A( :, : ) X REAL(WP), INTENT( OUT ) :: R(:), C(:) X REAL(WP), INTENT( OUT ), OPTIONAL :: ROWCND, COLCND, AMAX X END SUBROUTINE SGEEQU_F90 X SUBROUTINE DGEEQU_F90( A, R, C, ROWCND, COLCND, AMAX, INFO ) X USE LA_PRECISION, ONLY: WP => DP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X REAL(WP), INTENT( IN ) :: A( :, : ) X REAL(WP), INTENT( OUT ) :: R(:), C(:) X REAL(WP), INTENT( OUT ), OPTIONAL :: ROWCND, COLCND, AMAX X END SUBROUTINE DGEEQU_F90 X SUBROUTINE CGEEQU_F90( A, R, C, ROWCND, COLCND, AMAX, INFO ) X USE LA_PRECISION, ONLY: WP => SP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X COMPLEX(WP), INTENT( IN ) :: A( :, : ) X REAL(WP), INTENT( OUT ) :: R(:), C(:) X REAL(WP), INTENT( OUT ), OPTIONAL :: ROWCND, COLCND, AMAX X END SUBROUTINE CGEEQU_F90 X SUBROUTINE ZGEEQU_F90( A, R, C, ROWCND, COLCND, AMAX, INFO ) X USE LA_PRECISION, ONLY: WP => DP X INTEGER, INTENT(OUT), OPTIONAL :: INFO X COMPLEX(WP), INTENT( IN ) :: A( :, : ) X REAL(WP), INTENT( OUT ) :: R(:), C(:) X REAL(WP), INTENT( OUT ), OPTIONAL :: ROWCND, COLCND, AMAX X END SUBROUTINE ZGEEQU_F90 X END INTERFACE X! XEND MODULE LAPACK90_INTERFACES END_OF_FILE if test 45546 -ne `wc -c <'laf90mod.f90'`; then echo shar: \"'laf90mod.f90'\" unpacked with wrong size! fi # end of 'laf90mod.f90' fi if test -f 'pwcr_drv.f90' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'pwcr_drv.f90'\" else echo shar: Extracting \"'pwcr_drv.f90'\" \(4494 characters\) sed "s/^X//" >'pwcr_drv.f90' <<'END_OF_FILE' X!=========================================================================! X! file pwcr_drv.f90 ! X!=========================================================================! X! IMPROVED CYCLIC REDUCTION FOR SOLVING QUEUEING PROBLEMS ! X! by D.A. Bini and B. Meini ! X! (bini@dm.unipi.it meini@dm.unipi.it) ! X! Fortran 90 Program version 1.0, January 30, 1997 ! X! DRIVER PROGRAM ! X!=========================================================================! X! This is the driver program for solving the matrix equation ! X! A_1+X A_2+X^2 A_3+....X^(m-1) A_m-X=O (1) ! X! by means of pointwise cyclic reduction, where X is the (nb x nb) unknown! X! matrix and the (nb x nb) nonnegative matrices A_i, i=1,2,..., are such ! X! that A_1+A_2+...+A_m is column stochastic. ! X!=========================================================================! X! The program reads from the standard input the name of the file ! X! () containing the data. Then reads from the file ! X! the data and calls the subroutine PWCR solving the ! X! matrix equation (1). ! X! The computed solution is stored in the file ! X! it is written row-wise. ! X! The file must contain on each row the following data: ! X! 1-st row: NB ( integer ) = Block dimension ! X! 2-nd row: M ( integer ) = Number of blocks A_i ! X! 3-rd row: EPS ( real(kind(0.d0)) ) = Precision level ! X! The subsequent rows: L, I, J, VAL, where ! X! L : ( integer ) is the index of the block A_L ! X! I,J : ( integer ) are the indices of the (I,J)-th entry of the block ! X! A_L ! X! VAL : ( real(kind(0.d0)) ) is the value of the (I,J)-th entry of the ! X! block A_L ! X! The last row: 0,0,0,0.0d0 which denotes the end of the file. ! X! The values that are not reported in this file are intended to be zero. ! X!=========================================================================! XPROGRAM driver X USE pwcr_interface X IMPLICIT NONE X REAL(KIND(0.d0)), DIMENSION(:,:,:), POINTER :: a X REAL(KIND(0.d0)), DIMENSION(:,:), POINTER :: g X REAL(KIND(0.d0)) :: err, eps, val X INTEGER :: m, nb, i,j,l X CHARACTER(len=10) :: ifn X LOGICAL :: ex X WRITE(*,*)'-------------------------------------------------------' X WRITE(*,*)'IMPROVED CYCLIC REDUCTION FOR SOLVING QUEUEING PROBLEMS' X WRITE(*,*)' by D.A. Bini and B. Meini' X WRITE(*,*)' Fortran 90 Program version 1.0, January 30, 1997' X WRITE(*,*)'-------------------------------------------------------' X WRITE(*,*) X WRITE(*,'(A)', advance='NO')' Input-file name = ' X READ(*,*) ifn X INQUIRE(file=ifn, exist=ex) X IF(.NOT.ex)THEN X WRITE(*,*) X WRITE(*,*)'THE FILE ',ifn,'DOES NOT EXIST' X STOP X ENDIF X OPEN(unit=10,file=ifn,status='old') X J = INDEX(ifn,'.') X IF (J.EQ.0) J = INDEX(ifn,' ') X OPEN(unit=11,file=ifn(1:j-1)//'.out') X READ(10,*) nb X READ(10,*) m X READ(10,*) eps X ALLOCATE(a(nb,nb,m)) X ! reading data X a=0.0d0 X READ(10,*) l,i,j,val X reading : DO WHILE ( i/=0 ) X a(i,j,l)=val X READ(10,*)l,i,j,val X END DO reading X CALL pwcr(a,eps,g,err) X WRITE(11,*)'-------------------------------------------------------' X WRITE(11,*)'IMPROVED CYCLIC REDUCTION FOR SOLVING QUEUEING PROBLEMS' X WRITE(11,*)' by D.A. Bini and B. Meini' X WRITE(11,*)' Fortran 90 Program version 1.0, January 30, 1997' X WRITE(11,*)'-------------------------------------------------------' X WRITE(11,*) X WRITE(*,*)'Residual error=',err X WRITE(11,*)'Residual error=',err X WRITE(*,*) X WRITE(11,*) X loopw1 : DO i=1,nb X loopw2 : DO j=1,nb X WRITE(*,*)g(i,j) X WRITE(11,*)g(i,j) X END DO loopw2 X END DO loopw1 X STOP XEND PROGRAM driver X X X X X X X X X X END_OF_FILE if test 4494 -ne `wc -c <'pwcr_drv.f90'`; then echo shar: \"'pwcr_drv.f90'\" unpacked with wrong size! fi # end of 'pwcr_drv.f90' fi if test -f 'pwcr_fft.f90' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'pwcr_fft.f90'\" else echo shar: Extracting \"'pwcr_fft.f90'\" \(26485 characters\) sed "s/^X//" >'pwcr_fft.f90' <<'END_OF_FILE' X!=========================================================================! X! file pwcr_fft.f90 ! X!=========================================================================! X! IMPROVED CYCLIC REDUCTION FOR SOLVING QUEUEING PROBLEMS ! X! by D.A. Bini and B. Meini ! X! (bini@dm.unipi.it meini@dm.unipi.it) ! X! Fortran 90 Program version 1.0, January 30 1997 ! X! FFT subroutines ! X!=========================================================================! X! This file contains a set of subroutines needed in order ! X! to perform computations involving the Discrete Fourier Transform (DFT). ! X! More precisely the following subroutines are reported: ! X! ! X! fillroots : compute the roots of 1 ! X! ifft1 : compute the inverse DFT ! X! fft1 : compute the DFT ! X! iffts1 : compute the inverse DFT of the DFT of a real vector ! X! (Problem 2) ! X! ffts1 : compute the DFT of a real vector (Problem 1) ! X! iffts2 : solve Problem 4 ! X! ffts2 : solve Problem 3 ! X! twiddle : scale a complex vector ! X! itwiddle : scale a complex vector ! X! ftb1 : compute the block DFT of a real block vector (Problem 1 for ! X! matrix polynomials) ! X! ftb2 : solve Problem 3 for matrix polynomials ! X! iftb1 : compute the block inverse DFT of the DFT of a real block ! X! vector (Problem 2 for matrix polynomials) ! X! iftb2 : solve Problem 4 for matrix polynomials ! X!=========================================================================! X! SUBROUTINE FILLROOTS ! X!=========================================================================! X! Compute the n-th roots of 1 if they have not been yet computed before, ! X! otherwise return. ! X! The real and the imaginary parts of the roots are stored in the vectors ! X! wr, wi, respectively. These vectors are put in a common block. ! X!=========================================================================! Xsubroutine fillroots(n) X implicit none X real(kind(0.d0)), dimension(:), pointer :: wr, wi X real(kind(0.d0)), dimension(:), allocatable, save :: wwr, wwi X real(kind(0.d0)) :: pi, pi2 X integer :: i, j, k, m, mi1, n X common wr, wi X intrinsic size, allocated X pi= 6.28318530717958647692528676656d0/n X k=size(wr,1) X if(.not.allocated(wwr))then X allocate(wwr(n)) X allocate(wwi(n)) X allocate(wr(n)) X allocate(wi(n)) X loop1 : do i=1,n X pi2=(i-1)*pi X wwr(i)=cos(pi2) X wwi(i)=sin(pi2) X end do loop1 X wr=wwr X wi=wwi X return X end if X if(n<=k)then X return X end if X deallocate(wr) X deallocate(wi) X allocate(wr(n)) X allocate(wi(n)) X m=n/k X loop2 : do i=1,k X mi1=m*(i-1) X wr(mi1+1)=wwr(i) X wi(mi1+1)=wwi(i) X loop3 : do j=2,m X wr(mi1+j)=cos(pi*(mi1+j-1)) X wi(mi1+j)=sin(pi*(mi1+j-1)) X end do loop3 X end do loop2 X deallocate(wwr) X deallocate(wwi) X allocate(wwr(n)) X allocate(wwi(n)) X wwr=wr X wwi=wi X return Xend subroutine fillroots X X X X X!========================================================================! X! SUBROUTINE IFFT1 ! X!========================================================================! X! Compute the Inverse Discrete Fourier Transform of the vector having ! X! real parts stored in the vector x and imaginary parts stored in the ! X! vector y. On output x and y contain the real and the imaginary parts ! X! of the transformed vector, respectively. ! X! The algorithm is an adaptation of the split-radix FFT by ! X! Dhuamel-Hollman (Sorensen et al. IEEE trans. ! X! on Acoustics Speech and Signal Processing, ASSP-34, 1986). ! X!========================================================================! Xsubroutine ifft1(x,y) X use fft_interface, only : fillroots X implicit none X real(kind(0.d0)), dimension(:), pointer :: x,y X real(kind(0.d0)), dimension(:), pointer :: wr,wi X real(kind(0.d0)) :: r1, r2, s1, s2, s3, ss1, & X ss3, cc1, cc3, xt, yt, un X integer :: i, j, k, m, n, nmax, n2, & X n4, ne, na, is, id, & X i0, i1, i2, i3, n1 ,na3 X common wr,wi X intrinsic size X n=size(x,1) X m=log(n*1.d0)/log(2.d0) X if(2**m=j)goto 101 X xt=x(j) X yt=y(j) X y(j)=y(i) X y(i)=yt X x(j)=x(i) X x(i)=xt X101 k=n/2 X102 if (k>=j)goto 103 X j=j-k X k=k/2 X goto 102 X103 j=j+k X end do loop14 X un=1.d0/n X x=x*un X y=y*un X return Xend subroutine ifft1 X X X X!========================================================================! X! SUBROUTINE FFT1 ! X!========================================================================! X! Compute the Discrete Fourier Transform of the vector having ! X! real parts stored in the vector x and imaginary parts stored in the ! X! vector y. On output x and y contain the real and the imaginary parts ! X! of the transformed vector, respectively. ! X! The algorithm is an adaptation of the split-radix FFT by ! X! Dhuamel-Hollman (Sorensen et al. IEEE trans. ! X! on Acoustics Speech and Signal Processing, ASSP-34, 1986). ! X!========================================================================! Xsubroutine fft1(x,y) X use fft_interface, only : fillroots X implicit none X real(kind(0.d0)), dimension(:), pointer :: x,y X real(kind(0.d0)), dimension(:), pointer :: wr,wi X real(kind(0.d0)) :: r1, r2, r3, s2, ss1, ss3, & X cc1, cc3, xt, yt, s1 X integer :: i, j, k, m, n, nmax, n1, & X n2, n4, ne, na, is, id, & X i0, i1, i2, i3, na3 X common wr,wi X intrinsic size X n=size(x,1) X m=log(n*1.d0)/log(2.d0) X if(2**m=j)goto 101 X xt=x(j) X yt=y(j) X y(j)=y(i) X y(i)=yt X x(j)=x(i) X x(i)=xt X101 k=n/2 X102 if (k>=j)goto 103 X j=j-k X k=k/2 X goto 102 X103 j=j+k X end do loop104 X ! ---------------LENGTH TWO BUTTERFLIES-------------- X is=1 X id=4 X70 loop60 : do i0=is,n,id X i1=i0+1 X r1=x(i0) X x(i0)=r1+x(i1) X x(i1)=r1-x(i1) X r1=-y(i0) X y(i0)=-r1+y(i1) X y(i1)=-r1-y(i1) X end do loop60 X is=2*id-1 X id=4*id X if (is'pwcr_int.f90' <<'END_OF_FILE' X!=========================================================================! X! file pwcr_int.f90 ! X!=========================================================================! X! IMPROVED CYCLIC REDUCTION FOR SOLVING QUEUEING PROBLEMS ! X! by D.A. Bini and B. Meini ! X! (bini@dm.unipi.it meini@dm.unipi.it) ! X! Fortran 90 Program version 1.0, January 30 1997 ! X! Interface File for PWCR subroutines ! X!=========================================================================! X! Interface file for the subroutines: ! X! pwcr ! X! computeG ! X! residual ! X! schur ! X! test ! X! sc1p ! X! sc2p ! X! scc2p ! X! sc1 ! X! sc2 ! X! prodc ! X! means ! X! pmeans ! X! scc2 ! X! solver ! X! solvec ! X!=========================================================================! XMODULE pwcr_interface X INTERFACE X SUBROUTINE pwcr(a,eps,g,err) X IMPLICIT NONE X REAL(KIND(0.d0)), DIMENSION(:,:,:), POINTER :: a X REAL(KIND(0.d0)), DIMENSION(:,:), POINTER :: g X REAL(KIND(0.d0)) :: err, eps X END SUBROUTINE pwcr X endinterface X !======================================================================= X INTERFACE X SUBROUTINE computeG(a0,a0s,ae,ao,hae,hao,eps,g) X IMPLICIT NONE X REAL(KIND(0.d0)), DIMENSION(:,:,:), POINTER :: ae, ao, hae, hao X REAL(KIND(0.d0)), DIMENSION(:,:), POINTER :: a0, g X REAL(KIND(0.d0)), DIMENSION(:), POINTER :: a0s X REAL(KIND(0.d0)) :: eps X END SUBROUTINE computeG X endinterface X !======================================================================= X INTERFACE X SUBROUTINE residual(a,g,err) X IMPLICIT NONE X REAL(KIND(0.d0)), DIMENSION(:,:,:), POINTER :: a X REAL(KIND(0.d0)), DIMENSION(:,:), POINTER :: g X REAL(KIND(0.d0)) :: err X END SUBROUTINE residual X endinterface XEND MODULE pwcr_interface X!======================================================================= X!======================================================================= XMODULE schur_interface X INTERFACE X SUBROUTINE schur(ae,ao,hae,hao,mean,hmean,a1s,eps) X IMPLICIT NONE X REAL(KIND(0.d0)), DIMENSION(:,:,:), POINTER :: ae, ao, hae, hao X REAL(KIND(0.d0)), DIMENSION(:), POINTER :: mean, hmean, a1s X REAL(KIND(0.d0)) :: eps X END SUBROUTINE schur X END INTERFACE X !======================================================================= X INTERFACE X SUBROUTINE test(u,mean,eps,answer) X IMPLICIT NONE X REAL(KIND(0.d0)), POINTER, DIMENSION(:,:,:) :: u X REAL(KIND(0.d0)), POINTER, DIMENSION(:) :: mean X REAL(KIND(0.d0)) :: eps X LOGICAL :: answer X END SUBROUTINE test X END INTERFACE X !======================================================================= X INTERFACE X SUBROUTINE sc1p(tae,tao,a) X IMPLICIT NONE X REAL(KIND(0.d0)), POINTER, DIMENSION(:,:,:) :: tao, tae, a X END SUBROUTINE sc1p X END INTERFACE X !======================================================================= X INTERFACE X SUBROUTINE sc2p(tae,tao,a1,a2) X IMPLICIT NONE X REAL(KIND(0.d0)), POINTER, DIMENSION(:,:,:) :: tao, tae, a1, a2 X END SUBROUTINE sc2p X END INTERFACE X !======================================================================= X INTERFACE X SUBROUTINE scc2p(thae,thao,a1,a2) X IMPLICIT NONE X REAL(KIND(0.d0)), POINTER, DIMENSION(:,:,:) :: thao, thae, a1, a2 X END SUBROUTINE scc2p X END INTERFACE X !======================================================================= X INTERFACE X SUBROUTINE sc1(tae,tao,a) X IMPLICIT NONE X REAL(KIND(0.d0)), POINTER, DIMENSION(:,:,:) :: tae, tao, a X END SUBROUTINE sc1 X END INTERFACE X !======================================================================= X INTERFACE X SUBROUTINE sc2(tae,tao,a1,a2) X IMPLICIT NONE X REAL(KIND(0.d0)), POINTER, DIMENSION(:,:,:) :: tae, tao, a1, a2 X END SUBROUTINE sc2 X END INTERFACE X !======================================================================= X INTERFACE X SUBROUTINE prodc(a,b,ir,ic,rmat,cmat) X IMPLICIT NONE X REAL(KIND(0.d0)), POINTER, DIMENSION(:,:,:) :: a, b X REAL(KIND(0.d0)), POINTER, DIMENSION(:,:) :: rmat, cmat X INTEGER :: ir,ic X END SUBROUTINE prodc X END INTERFACE X !======================================================================= X INTERFACE X SUBROUTINE means(tae,tao,thae,mean,hmean) X IMPLICIT NONE X REAL(KIND(0.d0)), POINTER, DIMENSION(:,:,:) :: tae,tao, thae X REAL(KIND(0.d0)), POINTER, DIMENSION(:) :: mean,hmean X END SUBROUTINE means X END INTERFACE X !======================================================================= X INTERFACE X SUBROUTINE pmeans(ae,ao,hae,hao,mean,hmean) X IMPLICIT NONE X REAL(KIND(0.d0)), POINTER, DIMENSION(:,:,:) :: ae,ao, hae, hao X REAL(KIND(0.d0)), POINTER, DIMENSION(:) :: mean,hmean X END SUBROUTINE pmeans X END INTERFACE X !======================================================================= X INTERFACE X SUBROUTINE scc2(thae,thao,a1,a2) X IMPLICIT NONE X REAL(KIND(0.d0)), POINTER, DIMENSION(:,:,:) :: thae,thao, a1, a2 X END SUBROUTINE scc2 X END INTERFACE X !======================================================================= X INTERFACE X SUBROUTINE solver(rmat,termr) X IMPLICIT NONE X REAL(KIND(0.d0)), POINTER, DIMENSION(:,:) :: rmat,termr X END SUBROUTINE solver X END INTERFACE X !======================================================================= X INTERFACE X SUBROUTINE solvec(cmat,termc) X IMPLICIT NONE X COMPLEX(KIND(0.d0)), POINTER, DIMENSION(:,:) :: cmat,termc X END SUBROUTINE solvec X END INTERFACE X XEND MODULE schur_interface X X X X X X X X X X X X X X END_OF_FILE if test 7230 -ne `wc -c <'pwcr_int.f90'`; then echo shar: \"'pwcr_int.f90'\" unpacked with wrong size! fi # end of 'pwcr_int.f90' fi if test -f 'pwcr_sub.f90' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'pwcr_sub.f90'\" else echo shar: Extracting \"'pwcr_sub.f90'\" \(45750 characters\) sed "s/^X//" >'pwcr_sub.f90' <<'END_OF_FILE' X!=========================================================================! X! file pwcr_sub.f90 ! X!=========================================================================! X! IMPROVED CYCLIC REDUCTION FOR SOLVING QUEUEING PROBLEMS ! X! by D.A. Bini and B. Meini ! X! (bini@dm.unipi.it meini@dm.unipi.it) ! X! Fortran 90 Program version 1.0, January 30, 1997 ! X! SUBROUTINES ! X!===============