C ALGORITHM 782, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 24,NO. 2, June, 1998, P. 254--257. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # rrqr_acm/ # rrqr_acm/README # rrqr_acm/cv15.acm/ # rrqr_acm/cv15.acm/GenCode # rrqr_acm/cv15.acm/Makefile.GenCode # rrqr_acm/cv15.acm/cgeqpb.f # rrqr_acm/cv15.acm/cgeqpc.f # rrqr_acm/cv15.acm/cgeqpw.f # rrqr_acm/cv15.acm/cgeqpx.f # rrqr_acm/cv15.acm/cgeqpy.f # rrqr_acm/cv15.acm/clasmx.f # rrqr_acm/cv15.acm/clauc1.f # rrqr_acm/cv15.acm/cmylap.f # rrqr_acm/cv15.acm/ctrqpx.f # rrqr_acm/cv15.acm/ctrqpy.f # rrqr_acm/cv15.acm/ctrqxc.f # rrqr_acm/cv15.acm/ctrqyc.f # rrqr_acm/cv15.acm/ctrrnk.f # rrqr_acm/cv15.acm/zgeqpb.f # rrqr_acm/cv15.acm/zgeqpc.f # rrqr_acm/cv15.acm/zgeqpw.f # rrqr_acm/cv15.acm/zgeqpx.f # rrqr_acm/cv15.acm/zgeqpy.f # rrqr_acm/cv15.acm/zlasmx.f # rrqr_acm/cv15.acm/zlauc1.f # rrqr_acm/cv15.acm/zmylap.f # rrqr_acm/cv15.acm/ztrqpx.f # rrqr_acm/cv15.acm/ztrqpy.f # rrqr_acm/cv15.acm/ztrqxc.f # rrqr_acm/cv15.acm/ztrqyc.f # rrqr_acm/cv15.acm/ztrrnk.f # rrqr_acm/generate # rrqr_acm/lib/ # rrqr_acm/lib/Makefile # rrqr_acm/matgen/ # rrqr_acm/matgen/Makefile # rrqr_acm/matgen/clagge.f # rrqr_acm/matgen/claghe.f # rrqr_acm/matgen/clagsy.f # rrqr_acm/matgen/clarge.f # rrqr_acm/matgen/clarnd.f # rrqr_acm/matgen/claror.f # rrqr_acm/matgen/clarot.f # rrqr_acm/matgen/clatm1.f # rrqr_acm/matgen/clatm2.f # rrqr_acm/matgen/clatm3.f # rrqr_acm/matgen/clatme.f # rrqr_acm/matgen/clatmr.f # rrqr_acm/matgen/clatms.f # rrqr_acm/matgen/dlagge.f # rrqr_acm/matgen/dlagsy.f # rrqr_acm/matgen/dlaran.f # rrqr_acm/matgen/dlarge.f # rrqr_acm/matgen/dlarnd.f # rrqr_acm/matgen/dlaror.f # rrqr_acm/matgen/dlarot.f # rrqr_acm/matgen/dlatm1.f # rrqr_acm/matgen/dlatm2.f # rrqr_acm/matgen/dlatm3.f # rrqr_acm/matgen/dlatme.f # rrqr_acm/matgen/dlatmr.f # rrqr_acm/matgen/dlatms.f # rrqr_acm/matgen/slagge.f # rrqr_acm/matgen/slagsy.f # rrqr_acm/matgen/slaran.f # rrqr_acm/matgen/slarge.f # rrqr_acm/matgen/slarnd.f # rrqr_acm/matgen/slaror.f # rrqr_acm/matgen/slarot.f # rrqr_acm/matgen/slatm1.f # rrqr_acm/matgen/slatm2.f # rrqr_acm/matgen/slatm3.f # rrqr_acm/matgen/slatme.f # rrqr_acm/matgen/slatmr.f # rrqr_acm/matgen/slatms.f # rrqr_acm/matgen/zlagge.f # rrqr_acm/matgen/zlaghe.f # rrqr_acm/matgen/zlagsy.f # rrqr_acm/matgen/zlarge.f # rrqr_acm/matgen/zlarnd.f # rrqr_acm/matgen/zlaror.f # rrqr_acm/matgen/zlarot.f # rrqr_acm/matgen/zlatm1.f # rrqr_acm/matgen/zlatm2.f # rrqr_acm/matgen/zlatm3.f # rrqr_acm/matgen/zlatme.f # rrqr_acm/matgen/zlatmr.f # rrqr_acm/matgen/zlatms.f # rrqr_acm/testing/ # rrqr_acm/testing/ctest.lg.in # rrqr_acm/testing/ctest.me.in # rrqr_acm/testing/ctest.sm.in # rrqr_acm/testing/dtest.lg.in # rrqr_acm/testing/dtest.me.in # rrqr_acm/testing/dtest.sm.in # rrqr_acm/testing/stest.lg.in # rrqr_acm/testing/stest.me.in # rrqr_acm/testing/stest.sm.in # rrqr_acm/testing/testall # rrqr_acm/testing/testall.lg # rrqr_acm/testing/testall.me # rrqr_acm/testing/testall.sm # rrqr_acm/testing/v2/ # rrqr_acm/testing/v2/GenCode # rrqr_acm/testing/v2/Makefile # rrqr_acm/testing/v2/Makefile.GenCode # rrqr_acm/testing/v2/aladhd.f # rrqr_acm/testing/v2/alaerh.f # rrqr_acm/testing/v2/alaesm.f # rrqr_acm/testing/v2/alahd.f # rrqr_acm/testing/v2/alareq.f # rrqr_acm/testing/v2/alasum.f # rrqr_acm/testing/v2/alasvm.f # rrqr_acm/testing/v2/cchkaa.f # rrqr_acm/testing/v2/cchkqp.f # rrqr_acm/testing/v2/cchkrr.f # rrqr_acm/testing/v2/cerrqp.f # rrqr_acm/testing/v2/cerrrr.f # rrqr_acm/testing/v2/chkxer.f # rrqr_acm/testing/v2/cqpt01.f # rrqr_acm/testing/v2/cqrt11.f # rrqr_acm/testing/v2/cqrt12.f # rrqr_acm/testing/v2/crrt01.f # rrqr_acm/testing/v2/crrt02.f # rrqr_acm/testing/v2/dchkaa.f # rrqr_acm/testing/v2/dchkqp.f # rrqr_acm/testing/v2/dchkrr.f # rrqr_acm/testing/v2/derrqp.f # rrqr_acm/testing/v2/derrrr.f # rrqr_acm/testing/v2/dlaord.f # rrqr_acm/testing/v2/dqpt01.f # rrqr_acm/testing/v2/dqrt11.f # rrqr_acm/testing/v2/dqrt12.f # rrqr_acm/testing/v2/drrt01.f # rrqr_acm/testing/v2/drrt02.f # rrqr_acm/testing/v2/ilaenv.f # rrqr_acm/testing/v2/schkaa.f # rrqr_acm/testing/v2/schkqp.f # rrqr_acm/testing/v2/schkrr.f # rrqr_acm/testing/v2/serrqp.f # rrqr_acm/testing/v2/serrrr.f # rrqr_acm/testing/v2/slaord.f # rrqr_acm/testing/v2/sqpt01.f # rrqr_acm/testing/v2/sqrt11.f # rrqr_acm/testing/v2/sqrt12.f # rrqr_acm/testing/v2/srrt01.f # rrqr_acm/testing/v2/srrt02.f # rrqr_acm/testing/v2/xerbla.f # rrqr_acm/testing/v2/xlaenv.f # rrqr_acm/testing/v2/zchkaa.f # rrqr_acm/testing/v2/zchkqp.f # rrqr_acm/testing/v2/zchkrr.f # rrqr_acm/testing/v2/zerrqp.f # rrqr_acm/testing/v2/zerrrr.f # rrqr_acm/testing/v2/zqpt01.f # rrqr_acm/testing/v2/zqrt11.f # rrqr_acm/testing/v2/zqrt12.f # rrqr_acm/testing/v2/zrrt01.f # rrqr_acm/testing/v2/zrrt02.f # rrqr_acm/testing/ztest.lg.in # rrqr_acm/testing/ztest.me.in # rrqr_acm/testing/ztest.sm.in # rrqr_acm/timing/ # rrqr_acm/timing/ctime.lg.in # rrqr_acm/timing/ctime.me.in # rrqr_acm/timing/ctime.sm.in # rrqr_acm/timing/dtime.lg.in # rrqr_acm/timing/dtime.me.in # rrqr_acm/timing/dtime.sm.in # rrqr_acm/timing/stime.lg.in # rrqr_acm/timing/stime.me.in # rrqr_acm/timing/stime.sm.in # rrqr_acm/timing/timeall # rrqr_acm/timing/v2/ # rrqr_acm/timing/v2/GenCode # rrqr_acm/timing/v2/Makefile # rrqr_acm/timing/v2/Makefile.GenCode # rrqr_acm/timing/v2/atimck.f # rrqr_acm/timing/v2/atimin.f # rrqr_acm/timing/v2/ctimaa.f # rrqr_acm/timing/v2/ctimmg.f # rrqr_acm/timing/v2/ctimmm.f # rrqr_acm/timing/v2/ctimmv.f # rrqr_acm/timing/v2/ctimqp.f # rrqr_acm/timing/v2/ctimqr.f # rrqr_acm/timing/v2/ctimrr.f # rrqr_acm/timing/v2/dmflop.f # rrqr_acm/timing/v2/dopbl2.f # rrqr_acm/timing/v2/dopbl3.f # rrqr_acm/timing/v2/dopla.f # rrqr_acm/timing/v2/dprtb4.f # rrqr_acm/timing/v2/dprtb5.f # rrqr_acm/timing/v2/dprtbl.f # rrqr_acm/timing/v2/dtimaa.f # rrqr_acm/timing/v2/dtimmg.f # rrqr_acm/timing/v2/dtimmm.f # rrqr_acm/timing/v2/dtimmv.f # rrqr_acm/timing/v2/dtimqp.f # rrqr_acm/timing/v2/dtimqr.f # rrqr_acm/timing/v2/dtimrr.f # rrqr_acm/timing/v2/icopy.f # rrqr_acm/timing/v2/ilaenv.f # rrqr_acm/timing/v2/smflop.f # rrqr_acm/timing/v2/sopbl2.f # rrqr_acm/timing/v2/sopbl3.f # rrqr_acm/timing/v2/sopla.f # rrqr_acm/timing/v2/sprtb4.f # rrqr_acm/timing/v2/sprtb5.f # rrqr_acm/timing/v2/sprtbl.f # rrqr_acm/timing/v2/stimaa.f # rrqr_acm/timing/v2/stimmg.f # rrqr_acm/timing/v2/stimmm.f # rrqr_acm/timing/v2/stimmv.f # rrqr_acm/timing/v2/stimqp.f # rrqr_acm/timing/v2/stimqr.f # rrqr_acm/timing/v2/stimrr.f # rrqr_acm/timing/v2/xlaenv.f # rrqr_acm/timing/v2/ztimaa.f # rrqr_acm/timing/v2/ztimmg.f # rrqr_acm/timing/v2/ztimmm.f # rrqr_acm/timing/v2/ztimmv.f # rrqr_acm/timing/v2/ztimqp.f # rrqr_acm/timing/v2/ztimqr.f # rrqr_acm/timing/v2/ztimrr.f # rrqr_acm/timing/ztime.lg.in # rrqr_acm/timing/ztime.me.in # rrqr_acm/timing/ztime.sm.in # rrqr_acm/v15.acm/ # rrqr_acm/v15.acm/Dqr.in # rrqr_acm/v15.acm/GenCode # rrqr_acm/v15.acm/Makefile # rrqr_acm/v15.acm/Makefile.GenCode # rrqr_acm/v15.acm/README # rrqr_acm/v15.acm/REVISIONS # rrqr_acm/v15.acm/Sqr.in # rrqr_acm/v15.acm/c # rrqr_acm/v15.acm/compile # rrqr_acm/v15.acm/dgeqpb.f # rrqr_acm/v15.acm/dgeqpc.f # rrqr_acm/v15.acm/dgeqpw.f # rrqr_acm/v15.acm/dgeqpx.f # rrqr_acm/v15.acm/dgeqpy.f # rrqr_acm/v15.acm/dgntst.f # rrqr_acm/v15.acm/dlasmx.f # rrqr_acm/v15.acm/dlauc1.f # rrqr_acm/v15.acm/dmylap.f # rrqr_acm/v15.acm/dqr.f # rrqr_acm/v15.acm/dqrmtx.f # rrqr_acm/v15.acm/dtrqpx.f # rrqr_acm/v15.acm/dtrqpy.f # rrqr_acm/v15.acm/dtrqxc.f # rrqr_acm/v15.acm/dtrqyc.f # rrqr_acm/v15.acm/dtrrnk.f # rrqr_acm/v15.acm/dutils.f # rrqr_acm/v15.acm/esm.f # rrqr_acm/v15.acm/ilaenv.f # rrqr_acm/v15.acm/sgeqpb.f # rrqr_acm/v15.acm/sgeqpc.f # rrqr_acm/v15.acm/sgeqpw.f # rrqr_acm/v15.acm/sgeqpx.f # rrqr_acm/v15.acm/sgeqpy.f # rrqr_acm/v15.acm/sgntst.f # rrqr_acm/v15.acm/slasmx.f # rrqr_acm/v15.acm/slauc1.f # rrqr_acm/v15.acm/smylap.f # rrqr_acm/v15.acm/sqr.f # rrqr_acm/v15.acm/sqrmtx.f # rrqr_acm/v15.acm/strqpx.f # rrqr_acm/v15.acm/strqpy.f # rrqr_acm/v15.acm/strqxc.f # rrqr_acm/v15.acm/strqyc.f # rrqr_acm/v15.acm/strrnk.f # rrqr_acm/v15.acm/sutils.f # This archive created: Tue Jan 19 19:19:54 1999 export PATH; PATH=/bin:$PATH if test ! -d 'rrqr_acm' then mkdir 'rrqr_acm' fi cd 'rrqr_acm' if test -f 'README' then echo shar: will not over-write existing file "'README'" else cat << SHAR_EOF > 'README' ******************************************************************************* * * * RRQR Factorization * * "Codes for Rank-Revealing Factorizations of Dense Matrices" * * by C.H.Bischof and G.Quintana-Orti. * * * ******************************************************************************* * RELEASE: 1.0.5 (5-Jan-1997) * ******************************************************************************* Requirements: ============= To use these codes, both LAPACK and BLAS libraries are strictly necessary and must be correctly installed in the target machine. In addition, tuned BLAS libraries are strongly recommended to obtain optimal performance. Contents of the distribution file: ================================== The distribution file, named "rrqr_acm.tar.gz", contains all the directories and files needed for creating, testing and timings the new RRQR codes. The untaring of the file "rrqr_acm.tar.gz" creates the directory "rrqr_acm" in the current directory. All the other directories and files needed are created inside "rrqr_acm". So this file must be moved to the final position before installing it. Two steps are needed to uncompress and untar the code: 1. Uncompressing the file: gunzip rrqr_acm.tar.gz 2. Untaring the files and directories: tar xf rrqr_acm.tar Once uncompressed and untared, the directory "rrqr_acm" appears in the current directory. It will contain 6 subdirectories: lib: This directory contain the tool to create the RRQR object library. testing: This directory contains the drivers to test REAL, DOUBLE PRECISION, COMPLEX and COMPLEX*16 RRQR code. timing: This directory contains the drivers to get the timings for REAL, DOUBLE PRECISION, COMPLEX and COMPLEX*16 RRQR code. v15.acm: This directory contains the latest version of the REAL RRQR code. Besides, it contains a complete testing/timing driver for the RRQR code. This is the one used to get the timings shown in the paper. cv15.acm: This directory contains the latest version of the COMPLEX RRQR code. matgen: This directory is only needed if library "tmglib.a" is not installed. The library "tmglib.a" (timing library) is the timing library that comes with the LAPACK installation package. It contains some interesting matrix generators needed by both our testing and timing codes. If the library "tmglib.a" is already installed in your target machine, you can remove this directory. If you do not have this library, you can build it from this directory. To continue the installation, jump to "rrqr_acm" directory: cd rrqr_acm Steps to build the RRQR library: ================================ 1. Go to "lib" directory: cd lib 2. Edit the file "Makefile" and set the two macros FORTRAN and OPTS according to your system requirements. You can find this macros at the very beginning of the file. For example, on SUN the macros should be set to: FORTRAN = f77 # Fortran-77 compiler command. OPTS = -O # Fortran-77 compiler optimizing options. You might need to set other macros such as RANLIB, AR, etc. Ask your system administrator about these ones (usually well known in UNIX environments). 3. Type "make" to build the library for the four data-types: single, double, complex and complex*16. If not all the data-types are needed, type: "make single" to build the library for only REAL code, "make double" to build the library for only DOUBLE PRECISION code, "make complex" to build the library for only COMPLEX code, or "make complex*16" to build the library for only COMPLEX*16 code, or whatever combination you need, such as "make single complex", etc. This step will build the object files and the library. The object files are left in the directory "rrqr_acm/lib". The library is left in the directory "rrqr_acm" (the parent of directory "lib"), so as to allow an easy access from the testing and timing drivers. Steps to build the testing drivers: =================================== 1. Go to the testing directory: cd testing/v2 2. Edit the file "Makefile" and set the macros accordingly to your system. The macros you must set are the following: FORTRAN, OPTS, LOADER, LOADOPTS, LAPACKLIB, TMGLIB, and BLASLIB. You can find them at the very beginning of this file. For example, on SUN the macros should be set to: FORTRAN = f77 # Fortran-77 compiler command. OPTS = -O # Fortran-77 compiler optimizing options. LOADER = f77 # Fortran-77 linker command. LOADOPTS = -O # Fortran-77 linker optimizing options. LAPACKLIB = /usr/lib/lapack.a # Path of LAPACK library. TMGLIB = /usr/lib/tmglib.a # Path of TMGLIB library. BLASLIB = /usr/lib/blas.a # Path of BLAS library. You may need to ask your system administrator the paths of lapack, tmglib, and blas libraries. 3. The command "make" will make the testing programs for the four data-types: REAL, DOUBLE PRECISION, COMPLEX and COMPLEX*16. If not all the data-types are needed, type: "make single" to build the library for only REAL code, "make double" to build the library for only DOUBLE PRECISION code, "make complex" to build the library for only COMPLEX code, or "make complex*16" to build the library for only COMPLEX*16 code, or whatever combination you need, such as "make single complex", etc. The executable files are created in the directory "testing". 4. Now you can start the tests. If you want to perform some complete tests on you machine, just type "testall". This script performs some thorough and long tests for the four data types on three sets of matrices: small (files _test.sm.in), medium (files _test.me.in), and large (files _test.lg.in). In this case, the results will be placed in files _test.sm.out, _test.me.out and _test.lg.out. To perform some specific tests, you must type: xlintst_ < _test.in where _ must be replaced by s, d, c, or z, according to the data types you want to test. The input files (such as _test.in) tells the driver which matrix sizes, matrix blocks, matrix types, etc. to test. You can find many examples of input files in directory "testing": "stest.sm.in" is an input file to perform some small tests on single precision real numbers. "dtest.me.in" is an input file to perform some medium tests on double precision real numbers. "ctest.lg.in" is an input file to perform some large tests on single precision complex numbers.g etc. Note: Be careful with matrix sizes. If you are going to test big matrices, you must set the parameters of files "xchkaa.f" accordingly. The current drivers are prepared to test up to 500x500 matrices. If you want to test larger ones, you must change the definition of NMAX in the mentioned files. Steps to build the timing drivers: ================================== 1. Go to the timing directory: "cd timing/v2" 2. Edit the file "Makefile" and set the macros accordingly to your system. The macros you must set are the following: FORTRAN, OPTS, LOADER, LOADOPTS, LAPACKLIB, TMGLIB, and BLASLIB. You can find them at the very beginning of this file. For example, on SUN the macros should be set to: FORTRAN = f77 # Fortran-77 compiler command. OPTS = -O # Fortran-77 compiler optimizing options. LOADER = f77 # Fortran-77 linker command. LOADOPTS = -O # Fortran-77 linker optimizing options. LAPACKLIB = /usr/lib/lapack.a # Path of LAPACK library. TMGLIB = /usr/lib/tmglib.a # Path of TMGLIB library. BLASLIB = /usr/lib/blas.a # Path of BLAS library. You may need to ask your system administrator the paths of lapack, tmglib, and blas libraries. 3. The command "make" will make the timing programs for the four data-types. If not all the data-types are needed, type: "make single" to build the library for only REAL code, "make double" to build the library for only DOUBLE PRECISION code, "make complex" to build the library for only COMPLEX code, or "make complex*16" to build the library for only COMPLEX*16 code, or whatever combination you need, such as "make single complex", etc. The executable files are created in "timing" directory. 4. Now you can start the timings. If you want to perform some complete timings on you machine, just type "timeall". This script performs some thorough and long timings for the four data types on three sets of matrices: small (files _time.sm.in), medium (files _time.me.in), and large (files _time.lg.in). In this case, the results will be placed in files _time.sm.out, _time.me.out and _time.lg.out. To perform some specific timings, you must type: xlintim_ < _time.in where _ must be replaced by s, d, c, or z, according to the data types you want to test. The input files (such as _time.in) tells the driver which matrix sizes, matrix blocks, matrix types, etc. to test. You can find many examples of input files in directory "timing": "stime.sm.in" is an input file to perform some small timings on single precision real numbers. "dtime.me.in" is an input file to perform some medium timings on double precision real numbers. "ctime.lg.in" is an input file to perform some large timings on single precision complex numbers.g etc. Note: Be careful with matrix sizes. If you are going to test big matrices, you must set the parameters of files "xtimaa.f" accordingly. The current drivers are prepared to test up to 1001x1001 matrices. If you want to test larger ones, you must change the definition of NMAX in the mentioned files. Drivers in v15.acm: =================== 1. Go to the directory "v15.acm": "cd v15.acm" 2. Edit the file "Makefile" and set the macros accordingly to your system. The macros you must set are the following: P, FORTRAN, OPTS, LOADER, LOADOPTS, LAPACKLIB, TMGLIB, and BLASLIB. You can find them at the very beginning of this file. For example, on SUN the macros should be set to: P = s # Precision (s or d). FORTRAN = f77 # Fortran-77 compiler command. OPTS = -O # Fortran-77 compiler optimizing options. LOADER = f77 # Fortran-77 linker command. LOADOPTS = -O # Fortran-77 linker optimizing options. LAPACKLIB = /usr/lib/lapack.a # Path of LAPACK library. TMGLIB = /usr/lib/tmglib.a # Path of TMGLIB library. BLASLIB = /usr/lib/blas.a # Path of BLAS library. You may need to ask your system administrator the paths of lapack, tmglib, and blas libraries. 3. The command "make" will make the timing programs for data type specified in the makefile (s for single precision and d for double precision). The name of the executables are "sqr" for single precision and "dqr" for double precision. The executable files are created in the same directory. 4. Now you can start the timings. To perform some specific timings, you must type: sqr (for single precision), and dqr (for double precision). File "sqr" is controlled by the input file "Sqr.in" and file "dqr" is controlled by the input file "Dqr.in". These two "*.in" files specify the matrix sizes, block sizes, minimum time to benchmark, and other information to use during the execution. You can find examples of input files in directory "v15.acm". In this case, the name of the output file is composed of "time" plus a suffix named in the first line of the input file. Note: Be careful with matrix sizes. If you are going to test big matrices, you must set the parameters of files "xtimaa.f" accordingly. The current drivers are prepared to test up to 1001x1001 matrices. If you want to test larger ones, you must change the definition of NMAX in the mentioned files. Steps to build the matrix generator: ==================================== 1. Go to the matrix generating directory: "cd matgen" 2. Edit the file "Makefile" and set the macros FORTRAN and OPTS according to your system requirements. For example, on SUN the macros should be set to: FORTRAN = f77 # Fortran-77 compiler command. OPTS = -O # Fortran-77 compiler optimizing options. You might need to set other macros such as RANLIB, AR, etc. Ask your system administrator about these ones (usually well known in UNIX environments). 3. The command "make" will make the timing programs for the four data-types. If not all the data-types are needed, type: "make single" to build the library for only REAL code, "make double" to build the library for only DOUBLE PRECISION code, "make complex" to build the library for only COMPLEX code, or "make complex*16" to build the library for only COMPLEX*16 code, or whatever combination you need, such as "make single complex", etc. Known problems: =============== We have detected that IBM ESSL library can produce residuals a bit larger than usual when testing the code on 500x500 matrices. We have also checked that the residuals go back under the usual threshold when using the BLAS library obtained from Netlib, which can be obtained with LAPACK distribution. Subroutine xBDSQR from LAPACK 2.0 seems to have an error. The same subroutine from LAPACK 1.0b is working fine. This subroutine computes the singular value decomposition of a bidiagonal matrix. This bug have affected our drivers in directory "v15.acm" by showing in a few cases big residuals in the comparison of the singular values of A and R. In case of problems: ==================== You can send e-mail to: gquintan@dsic.upv.es gquintan@inf.uji.es SHAR_EOF fi # end of overwriting check if test ! -d 'cv15.acm' then mkdir 'cv15.acm' fi cd 'cv15.acm' if test -f 'GenCode' then echo shar: will not over-write existing file "'GenCode'" else cat << SHAR_EOF > 'GenCode' make source -f Makefile.GenCode SHAR_EOF fi # end of overwriting check if test -f 'Makefile.GenCode' then echo shar: will not over-write existing file "'Makefile.GenCode'" else cat << SHAR_EOF > 'Makefile.GenCode' ####################################################################### # # Makefile for generating COMPLEX single and double precision source. # # Authors: C.H.Bischof and G.Quintana-Orti # ####################################################################### ####################################################################### # # The user must set this options: # REAL_SOURCES = ../../v15 COMPLEX_SOURCES = ../../cv15 GENERATE = ../generate CPP = /lib/cpp CPPFLAGS = "-I$(REAL_SOURCES)" # ####################################################################### ####################################################################### # # Modules for Rank-Revealing QR: # C_RRQR_MODULES = \ cgeqpb.f cgeqpw.f cgeqpc.f \ cgeqpx.f ctrqpx.f ctrqxc.f \ cgeqpy.f ctrqpy.f ctrqyc.f \ ctrrnk.f clauc1.f clasmx.f \ cmylap.f Z_RRQR_MODULES = \ zgeqpb.f zgeqpw.f zgeqpc.f \ zgeqpx.f ztrqpx.f ztrqxc.f \ zgeqpy.f ztrqpy.f ztrqyc.f \ ztrrnk.f zlauc1.f zlasmx.f \ zmylap.f # # ####################################################################### source: complex complex16 complex: $(C_RRQR_MODULES) complex16: $(Z_RRQR_MODULES) clean: - rm -f *.f # # Rules for the source modules of RRQR. # cgeqpb.f: $(COMPLEX_SOURCES)/ygeqpb.F $(GENERATE) c $(COMPLEX_SOURCES)/ygeqpb.F $@ $(CPP) $(CPPFLAGS) cgeqpw.f: $(COMPLEX_SOURCES)/ygeqpw.F $(GENERATE) c $(COMPLEX_SOURCES)/ygeqpw.F $@ $(CPP) $(CPPFLAGS) cgeqpc.f: $(COMPLEX_SOURCES)/ygeqpc.F $(GENERATE) c $(COMPLEX_SOURCES)/ygeqpc.F $@ $(CPP) $(CPPFLAGS) cgeqpx.f: $(COMPLEX_SOURCES)/ygeqpx.F $(GENERATE) c $(COMPLEX_SOURCES)/ygeqpx.F $@ $(CPP) $(CPPFLAGS) ctrqpx.f: $(COMPLEX_SOURCES)/ytrqpx.F $(GENERATE) c $(COMPLEX_SOURCES)/ytrqpx.F $@ $(CPP) $(CPPFLAGS) ctrqxc.f: $(COMPLEX_SOURCES)/ytrqxc.F $(GENERATE) c $(COMPLEX_SOURCES)/ytrqxc.F $@ $(CPP) $(CPPFLAGS) cgeqpy.f: $(COMPLEX_SOURCES)/ygeqpy.F $(GENERATE) c $(COMPLEX_SOURCES)/ygeqpy.F $@ $(CPP) $(CPPFLAGS) ctrqpy.f: $(COMPLEX_SOURCES)/ytrqpy.F $(GENERATE) c $(COMPLEX_SOURCES)/ytrqpy.F $@ $(CPP) $(CPPFLAGS) ctrqyc.f: $(COMPLEX_SOURCES)/ytrqyc.F $(GENERATE) c $(COMPLEX_SOURCES)/ytrqyc.F $@ $(CPP) $(CPPFLAGS) ctrrnk.f: $(COMPLEX_SOURCES)/ytrrnk.F $(GENERATE) c $(COMPLEX_SOURCES)/ytrrnk.F $@ $(CPP) $(CPPFLAGS) clauc1.f: $(COMPLEX_SOURCES)/ylauc1.F $(GENERATE) c $(COMPLEX_SOURCES)/ylauc1.F $@ $(CPP) $(CPPFLAGS) clasmx.f: $(COMPLEX_SOURCES)/ylasmx.F $(GENERATE) c $(COMPLEX_SOURCES)/ylasmx.F $@ $(CPP) $(CPPFLAGS) cmylap.f: $(COMPLEX_SOURCES)/cmylap.f cp $(COMPLEX_SOURCES)/cmylap.f . zgeqpb.f: $(COMPLEX_SOURCES)/ygeqpb.F $(GENERATE) z $(COMPLEX_SOURCES)/ygeqpb.F $@ $(CPP) $(CPPFLAGS) zgeqpw.f: $(COMPLEX_SOURCES)/ygeqpw.F $(GENERATE) z $(COMPLEX_SOURCES)/ygeqpw.F $@ $(CPP) $(CPPFLAGS) zgeqpc.f: $(COMPLEX_SOURCES)/ygeqpc.F $(GENERATE) z $(COMPLEX_SOURCES)/ygeqpc.F $@ $(CPP) $(CPPFLAGS) zgeqpx.f: $(COMPLEX_SOURCES)/ygeqpx.F $(GENERATE) z $(COMPLEX_SOURCES)/ygeqpx.F $@ $(CPP) $(CPPFLAGS) ztrqpx.f: $(COMPLEX_SOURCES)/ytrqpx.F $(GENERATE) z $(COMPLEX_SOURCES)/ytrqpx.F $@ $(CPP) $(CPPFLAGS) ztrqxc.f: $(COMPLEX_SOURCES)/ytrqxc.F $(GENERATE) z $(COMPLEX_SOURCES)/ytrqxc.F $@ $(CPP) $(CPPFLAGS) zgeqpy.f: $(COMPLEX_SOURCES)/ygeqpy.F $(GENERATE) z $(COMPLEX_SOURCES)/ygeqpy.F $@ $(CPP) $(CPPFLAGS) ztrqpy.f: $(COMPLEX_SOURCES)/ytrqpy.F $(GENERATE) z $(COMPLEX_SOURCES)/ytrqpy.F $@ $(CPP) $(CPPFLAGS) ztrqyc.f: $(COMPLEX_SOURCES)/ytrqyc.F $(GENERATE) z $(COMPLEX_SOURCES)/ytrqyc.F $@ $(CPP) $(CPPFLAGS) ztrrnk.f: $(COMPLEX_SOURCES)/ytrrnk.F $(GENERATE) z $(COMPLEX_SOURCES)/ytrrnk.F $@ $(CPP) $(CPPFLAGS) zlauc1.f: $(COMPLEX_SOURCES)/ylauc1.F $(GENERATE) z $(COMPLEX_SOURCES)/ylauc1.F $@ $(CPP) $(CPPFLAGS) zlasmx.f: $(COMPLEX_SOURCES)/ylasmx.F $(GENERATE) z $(COMPLEX_SOURCES)/ylasmx.F $@ $(CPP) $(CPPFLAGS) zmylap.f: $(COMPLEX_SOURCES)/zmylap.f cp $(COMPLEX_SOURCES)/zmylap.f . SHAR_EOF fi # end of overwriting check if test -f 'cgeqpb.f' then echo shar: will not over-write existing file "'cgeqpb.f'" else cat << SHAR_EOF > 'cgeqpb.f' SUBROUTINE CGEQPB( JOB, M, N, K, A, LDA, C, LDC, JPVT, $ IRCOND, ORCOND, RANK, SVLUES, WORK, LWORK, $ RWORK, INFO ) * * This code is part of a release of the package for computing * rank-revealing QR Factorizations written by: * ================================================================== * Christian H. Bischof and Gregorio Quintana-Orti * Math. and Comp. Sci. Div. Departamento de Informatica * Argonne National Lab. Universidad Jaime I * Argonne, IL 60439 Campus P. Roja, 12071 Castellon * USA Spain * bischof@mcs.anl.gov gquintan@inf.uji.es * ================================================================== * $Revision: 1.42 $ * $Date: 96/12/30 16:59:36 $ * * .. Scalar Arguments .. INTEGER JOB, M, N, K, LDA, LDC, RANK, LWORK, INFO REAL IRCOND, ORCOND * .. * .. Array Arguments .. INTEGER JPVT( * ) COMPLEX A( LDA, * ), C( LDC, * ), WORK( * ) REAL SVLUES( 4 ), RWORK( * ) * .. * * Purpose * ======= * * CGEQPB computes a QR factorization * A*P = Q*[ R11 R12 ] * [ 0 R22 ] * of a real m by n matrix A. The permutation P is * chosen with the goal to reveal the rank of A by a * suitably dimensioned trailing submatrix R22 with norm(R22) * being small. * * Arguments * ========= * * JOB (input) INTEGER * The job to do: * = 1: The transformations needed in the * triangularization are only applied to matrix A. * Thus, matrix C is not updated. * = 2: The same transformations needed in the * triangularization of matrix A are applied to * matrix C from the left. * That is, if Q'*A*P=R, then C := Q'*C. * In this case, matrix C is m-by-k. * = 3: The transpose of the transformations needed * in the triangularization of matrix A are applied * to matrix C from the right. * That is, if Q'*A*P=R, then C := C*Q. * In this case, matrix C is k-by-m. * In these three cases, the permutations are always saved * into vector JPVT. * * M (input) INTEGER * The number of rows of matrices A. M >= 0. * If JOB=2, M is the number of rows of matrix C. * If JOB=3, M is the number of columns of matrix C. * * N (input) INTEGER * The number of columns of matrix A. N >= 0. * * K (input) INTEGER * It defines the dimension of matrix C. K >= 0. * If JOB=2, K is the number of columns of matrix C. * If JOB=3, K is the number of rows of matrix C. * * A (input/output) COMPLEX array, dimension (LDA,N) * On entry, the m by n matrix A. * On exit, the upper triangle of the array contains the * min(m,n) by n upper trapezoidal matrix R; the lower triangle * array is filled with zeros. * * LDA (input) INTEGER * The leading dimension of array A. LDA >= MAX(1,M). * * C (input/output) COMPLEX array, dimension * ( LDC, K ) if JOB=2. * ( LDC, M ) if JOB=3. * If argument JOB asks, all the transformations * applied to matrix A are also applied to matrix C. * * LDC (input) INTEGER * The leading dimension of array C. * If JOB=2, then LDC >= MAX(1,M). * If JOB=3, then LDC >= MAX(1,K). * * JPVT (output) INTEGER array, dimension (N) * JPVT(I) = K <==> Column K of A has been permuted into * position I in AP. * JPVT(1:RANK) contains the indices of the columns considered * linearly independent. * JPVT(RANK+1:N) contains the indices of the columns considered * linearly dependent from the previous ones. * * IRCOND (input) FLOATING_DECLARE * 1/IRCOND specifies an upper bound on the condition number * of R11. If IRCOND == 0, IRCOND = machine precision is chosen * as default. IRCOND must be >= 0. * * ORCOND (output) FLOATING_DECLARE * 1/ORCOND is an estimate for the condition number of R11. * * RANK (output) INTEGER * RANK is an estimate for the numerical rank of A with respect * to the threshold 1/IRCOND in the sense that * RANK = arg_max(cond_no(R(1:r,1:r))<1/IRCOND) * This may be an underestimate of the rank if the leading * columns were not well-conditioned. * * SVLUES (output) REAL array, dimension (4) * On exit, SVLUES contains estimates of some of the * singular values of the triangular factor R. * SVLUES(1): largest singular value of R(1:RANK,1:RANK) * SVLUES(2): smallest singular value of R(1:RANK,1:RANK) * SVLUES(3): smallest singular value of R(1:RANK+1,1:RANK+1) * SVLUES(4): smallest singular value of R * If the triangular factorization is a rank-revealing one * (which will be the case if the leading columns were well- * conditioned), then SVLUES(1) will also be an estimate * for the largest singular value of A, SVLUES(2) and SVLUES(3) * will be estimates for the RANK-th and (RANK+1)-st singular * value of A, and SVLUES(4) will be an estimate for the * smallest singular value of A. * By examining these values, one can confirm that the rank is * well defined with respect to the threshold chosen. * * WORK (workspace) COMPLEX array, dimension (LWORK) * On exit: WORK(1) is the size of the storage array needed * for optimal performance * * LWORK (input) INTEGER * The dimension of array WORK. * If JOB=1: * The unblocked strategy requires that: * LWORK >= 2*MN+N. * The block algorithm requires that: * LWORK >= 2*MN+N*NB. * If JOB<>1: * The unblocked strategy requires that: * LWORK >= 2*MN+MAX(K,N). * The block algorithm requires that: * LWORK >= 2*MN+NB*NB+NB*MAX(K,N). * Where MN = min(M,N) and NB is the maximum of blocksize * used within xGEQRF and blocksize used within xUNMQR. * In both cases, the minimum required workspace is the * one for the unblocked strategy. * * RWORK (workspace) REAL array, dimension (2*N) * * INFO (output) INTEGER * = 0: Successful exit * < 0: If INFO = -i, the i-th argument had an illegal value * * ===================================================================== * * .. Parameters .. REAL ZERO COMPLEX CZERO PARAMETER ( ZERO = 0.0E+0, CZERO = ( 0.0E+0, 0.0E+0 ) ) INTEGER INB, INBMIN, IXOVER PARAMETER ( INB = 1, INBMIN = 2, IXOVER = 3 ) * * Indices into the 'svlues' array. * INTEGER IMAX, IBEFOR, IAFTER, IMIN PARAMETER ( IMAX = 1, IBEFOR = 2, IAFTER = 3, IMIN = 4 ) * .. * .. Local Scalars .. INTEGER I, J, MN, ITEMP, KK, LACPTD, MVIDX, STREJ, $ ACCPTD, NB, LWSIZE, NLLITY, KB, WSIZE, WKMIN REAL SMIN, MXNM, RCOND LOGICAL BLOCK * .. * .. External Subroutines .. EXTERNAL XERBLA, CGEQPC, CGEQPW, $ CLARFT, CLARFB * .. * .. External Functions .. INTEGER ILAENV REAL SLAMCH, SLASMX EXTERNAL ILAENV, SLAMCH, SLASMX * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN, REAL, INT * .. * .. Executable Statements .. * MN = MIN( M, N ) * * Compute the minimum required complex workspace. * IF( JOB.EQ.1 ) THEN WKMIN = 2*MN + N ELSE WKMIN = 2*MN + MAX( N, K ) END IF * * Test input arguments * ==================== * INFO = 0 IF( ( JOB.LT.1 ).OR.( JOB.GT.3 ) ) THEN INFO = -1 ELSE IF( M.LT.0 ) THEN INFO = -2 ELSE IF( M.LT.0 ) THEN INFO = -3 ELSE IF( K.LT.0 ) THEN INFO = -4 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -6 ELSE IF( ( ( JOB.EQ.1 ).AND.( LDC.LT.1 ) ).OR. $ ( ( JOB.EQ.2 ).AND.( LDC.LT.MAX( 1, M ) ) ).OR. $ ( ( JOB.EQ.3 ).AND.( LDC.LT.MAX( 1, K ) ) ) ) THEN INFO = -8 ELSE IF( IRCOND.LT.ZERO ) THEN INFO = -10 ELSE IF( LWORK.LT.MAX( 1, WKMIN ) ) THEN INFO = -15 END IF IF( ( INFO .EQ. 0 .OR. INFO .EQ. -15 ).AND. LWORK.GE.1 ) THEN * * Compute the optimal complex workspace. * IF( JOB.EQ.1 ) THEN NB = ILAENV( INB, 'CGEQRF', ' ', M, N, 0, 0 ) WSIZE = 2*MN + MAX( 3*N, N*NB ) ELSE NB = MAX( ILAENV( INB, 'CGEQRF', ' ', M, N, 0, 0 ), $ ILAENV( INB, 'CUNMQR', ' ', M, N, 0, 0 ) ) WSIZE = MAX( 2*MN + MAX( N, K ), $ 2*MN + NB*NB + NB*MAX( N, K ) ) END IF WORK( 1 ) = REAL( WSIZE ) END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CGEQPB', -INFO ) RETURN END IF * * Initialization of vector JPVT. * DO 70 J = 1, N JPVT( J ) = J 70 CONTINUE * * Quick return if possible * IF( MN.EQ.0 ) THEN RANK = 0 ORCOND = ZERO SVLUES( IMAX ) = ZERO SVLUES( IBEFOR ) = ZERO SVLUES( IAFTER ) = ZERO SVLUES( IMIN ) = ZERO RETURN END IF * * Check whether Threshold for condition number was supplied. * If not, choose machine precision as default for RCOND. * IF( IRCOND.GT.ZERO ) THEN RCOND = IRCOND ELSE RCOND = SLAMCH( 'Epsilon' ) END IF * * Determine block size and whether to use blocked code at all * IF( LWORK.LT.WSIZE ) THEN IF( JOB.EQ.1 ) THEN NB = ( LWORK-2*MN )/N ELSE ITEMP = INT( SQRT( REAL( $ MAX( K, N )**2+4*LWORK-8*MN ) ) ) NB = ( ITEMP-MAX( K, N ) )/2 END IF END IF * BLOCK = ( ( NB.GT.1 ).AND. $ ( NB.GE.ILAENV( INBMIN, 'CGEQRF', ' ', M, N, 0, 0 ) ).AND. $ ( MN.GE.ILAENV( IXOVER, 'CGEQRF', ' ', M, N, 0, 0 ) ) ) * * The size of the pivot window is chosen to be NB + NLLITY * for the blocked algorithm. * NLLITY = MIN( MN, MAX( 10, NB/2+(N*5)/100 ) ) * * *************************************************** * * Move column with largest residual norm up front * * *************************************************** * CALL CGEQPC( JOB, M, N, K, A, LDA, C, LDC, 1, 0, $ RCOND, LACPTD, JPVT, WORK( 1 ), WORK( MN+1 ), $ SVLUES, MXNM, WORK( 2*MN+1 ), LWORK-2*MN, RWORK ) IF( LACPTD.EQ.1 ) THEN IF( LACPTD.EQ.MN ) THEN RANK = 1 ORCOND = SVLUES( IBEFOR )/SVLUES( IMAX ) GOTO 30 ELSE SMIN = SVLUES( IBEFOR ) END IF ELSE RANK = 0 ORCOND = ZERO SVLUES( IMAX ) = ZERO SVLUES( IBEFOR ) = ZERO SVLUES( IAFTER ) = ZERO SVLUES( IMIN ) = ZERO GOTO 30 END IF * * **************************** * * Factor remaining columns * * **************************** * IF( BLOCK ) THEN * * *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* * * Using blocked code with restricted pivoting strategy * * *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* * STREJ = N+1 KK = 2 * 10 IF( ( KK.GE.STREJ ).OR.( KK.GT.MN ) ) GOTO 20 * * invariant: A(:,KK) is the first column in currently * considered block column. * KB = MIN( NB, MIN( MN+1, STREJ )-KK ) * * The goal now is to find "KB" independent columns * among the remaining STREJ-KK not yet rejected columns. * LWSIZE = MIN( STREJ-KK, KB+NLLITY ) CALL CGEQPW( M, LWSIZE, KB, KK-1, LACPTD, A, LDA, JPVT, $ RCOND, WORK( MN+1 ), SMIN, MXNM, $ WORK( 1 ), WORK( 2*MN+1 ), RWORK ) IF( LACPTD.GT.0 ) THEN * * Accumulate Householder vectors in a block reflector. * CALL CLARFT( 'Forward', 'Columnwise', M-KK+1, $ LACPTD, A( KK, KK ), LDA, WORK( KK ), $ WORK( 2*MN+1 ), LACPTD ) * * Apply block reflector to A(KK:M,KK+LWSIZE:N). * IF( ( KK+LWSIZE ).LE.N ) THEN CALL CLARFB( 'Left', 'Conjugate Transpose', $ 'Forward', 'Columnwise', $ M-KK+1, N-KK-LWSIZE+1, LACPTD, $ A( KK, KK ), LDA, $ WORK( 2*MN+1 ), LACPTD, $ A( KK, KK+LWSIZE ), LDA, $ WORK( 2*MN+LACPTD*LACPTD+1 ), $ N-KK-LWSIZE+1 ) END IF * * Apply block reflector to the corresponding part * of matrix C. * IF( ( JOB.EQ.2 ).AND.( K.GT.0 ) ) THEN * * Apply it to matrix C(KK:M,1:K) from the left. * CALL CLARFB( 'Left', 'Conjugate Transpose', $ 'Forward', 'Columnwise', M-KK+1, K, $ LACPTD, A( KK, KK ), LDA, $ WORK( 2*MN+1 ), LACPTD, $ C( KK, 1 ), LDC, $ WORK( 2*MN+LACPTD*LACPTD+1 ), K ) ELSE IF( ( JOB.EQ.3 ).AND.( K.GT.0 ) ) THEN * * Apply the transpose of it to matrix C(1:K,KK:M) * from the right. * CALL CLARFB( 'Right', 'No Transpose', 'Forward', $ 'Columnwise', K, M-KK+1, LACPTD, $ A( KK, KK ), LDA, $ WORK( 2*MN+1 ), LACPTD, $ C( 1, KK ), LDC, $ WORK( 2*MN+LACPTD*LACPTD+1 ), K ) END IF END IF * * Move rejected columns to the end if there is space. * IF( LACPTD.LT.KB ) THEN IF( STREJ.LE.( KK+LWSIZE ) ) THEN STREJ = KK + LACPTD ELSE MVIDX = STREJ DO 40 I = KK+LACPTD, $ MIN( KK+LWSIZE-1, STREJ-LWSIZE+LACPTD-1 ) MVIDX = MVIDX - 1 CALL CSWAP( M, A( 1, I ),1, A( 1, MVIDX ),1 ) ITEMP = JPVT( I ) JPVT( I ) = JPVT( MVIDX ) JPVT( MVIDX ) = ITEMP 40 CONTINUE STREJ = MVIDX END IF END IF KK = KK + LACPTD GOTO 10 20 CONTINUE ACCPTD = KK-1 SVLUES( IMAX ) = SLASMX( ACCPTD )*MXNM SVLUES( IBEFOR ) = SMIN IF( ACCPTD.LT.MN ) THEN * * Process rejected columns. * CALL CGEQPC( JOB, M, N, K, A, LDA, C, LDC, MN-KK+1, $ KK-1, RCOND, LACPTD, JPVT, WORK( 1 ), $ WORK( MN+1 ), SVLUES, MXNM, WORK( 2*MN+1 ), $ LWORK-2*MN, RWORK ) RANK = ACCPTD + LACPTD ELSE RANK = ACCPTD SVLUES( IAFTER ) = SMIN SVLUES( IMIN ) = SMIN END IF ELSE * * *-*-*-*-*-*-*-*-*-*-*-*-* * * using unblocked code * * *-*-*-*-*-*-*-*-*-*-*-*-* * ACCPTD = 1 CALL CGEQPC( JOB, M, N, K, A, LDA, C, LDC, MN-ACCPTD, $ ACCPTD, RCOND, LACPTD, JPVT, WORK( 1 ), $ WORK( MN+1 ), SVLUES, MXNM, WORK( 2*MN+1 ), $ LWORK-2*MN, RWORK ) RANK = ACCPTD+LACPTD * END IF ORCOND = SVLUES( IBEFOR )/SVLUES( IMAX ) * * Nullify the lower part of matrix A. * 30 CONTINUE DO 50 J = 1, MN DO 60 I = J+1, M A( I, J ) = CZERO 60 CONTINUE 50 CONTINUE * WORK( 1 ) = REAL( WSIZE ) RETURN * * End of CGEQPB * END SHAR_EOF fi # end of overwriting check if test -f 'cgeqpc.f' then echo shar: will not over-write existing file "'cgeqpc.f'" else cat << SHAR_EOF > 'cgeqpc.f' SUBROUTINE CGEQPC( JOB, M, N, K, A, LDA, C, LDC, DSRD, OFFSET, $ IRCOND, LACPTD, JPVT, TAU, X, SVLUES, MXNM, $ WORK, LWORK, RWORK ) * * This code is part of a release of the package for computing * rank-revealing QR Factorizations written by: * ================================================================== * Christian H. Bischof and Gregorio Quintana-Orti * Math. and Comp. Sci. Div. Departamento de Informatica * Argonne National Lab. Universidad Jaime I * Argonne, IL 60439 Campus P. Roja, 12071 Castellon * USA Spain * bischof@mcs.anl.gov gquintan@inf.uji.es * ================================================================== * $Revision: 1.42 $ * $Date: 96/12/30 16:59:37 $ * * .. Scalar Arguments .. INTEGER JOB, M, N, K, LDA, LDC, DSRD, OFFSET, LACPTD, $ LWORK REAL IRCOND, MXNM * .. * .. Array Arguments .. INTEGER JPVT( * ) COMPLEX A( LDA, * ), C( LDC, * ), TAU( * ), $ WORK( * ), X( * ) REAL SVLUES( 4 ), RWORK( * ) * .. * * Purpose: * ======= * * CGEQPC continues a partial QR factorization of A. If * A(1:OFFSET,1:OFFSET) has been reduced to upper triangular * form, then SGQPC applies the traditional column pivoting * strategy to identify DSRD more independent columns of A with * the restriction that the condition number of the leading * triangle of A should not be larger than 1/IRCOND. If * LACPTD ( <= DSRD) such columns are found, then the condition * number of * A(1:OFFSET+LACPTD,1:OFFSET+LACPTD) is less than 1/IRCOND. * If LACPTD < DSRD, then the QR factorization of A is completed, * otherwise only DSRD new steps were performed. * * Arguments: * ========= * * JOB (input) INTEGER * The job to do: * = 1: The transformations needed in the * triangularization are only applied to matrix A. * Thus, matrix C is not updated. * = 2: The same transformations needed in the * triangularization of matrix A are applied to * matrix C from the left. * That is, if Q'*A*P=R, then C := Q'*C. * In this case, matrix C is m-by-k. * = 3: The transpose of the transformations needed * in the triangularization of matrix A are applied * to matrix C from the right. * That is, if Q'*A*P=R, then C := C*Q. * In this case, matrix C is k-by-m. * In these three cases, the permutations are always saved * into vector JPVT. * * M (input) INTEGER * The number of rows of matrices A. M >= 0. * If JOB=2, M is the number of rows of matrix C. * If JOB=3, M is the number of columns of matrix C. * * N (input) INTEGER * The number of columns of matrix A. N >= 0. * * K (input) INTEGER * It defines the dimension of matrix C. K >= 0. * If JOB=2, K is the number of columns of matrix C. * If JOB=3, K is the number of rows of matrix C. * * A (input/output) COMPLEX array, dimension (LDA,N) * On entry, the m by n matrix A. * On exit, the upper triangle of the array contains the * min(m,n) by n upper trapezoidal matrix R; the lower triangle * array is filled with zeros. * * LDA (input) INTEGER * The leading dimension of array A. LDA >= max(1,M). * * C (input/output) COMPLEX array, dimension * ( LDC, K ) if JOB=2. * ( LDC, M ) if JOB=3. * If argument JOB asks, all the transformations * applied to matrix A are also applied to matrix C. * * LDC (input) INTEGER * The leading dimension of array C. * If JOB=2, then LDC >= MAX(1,M). * If JOB=3, then LDC >= MAX(1,K). * * DSRD (input) INTEGER * The number of independent columns one would like to * extract. * * OFFSET (input) INTEGER * A(1:OFFSET,1:OFFSET) has already been factored. * OFFSET >= 0. * * IRCOND (input) REAL * 1/IRCOND is threshold for condition number. * * LACPTD (output) INTEGER * The number of additional columns that were identified * as independent. * * JPVT (input/output) INTEGER array, dimension (N) * If JPVT(I) = K, then the Ith column of the permuted * A was the Kth column of the original A. * * TAU (input/output) COMPLEX array, dimension (MIN(M,N)) * Further details of the matrix Q (see A). * * X (input/output) COMPLEX array, dimension (MIN(M,N)) * On entry: X(1:OFFSET) contains an approximate smallest * left singular vector of A(1:OFFSET,1:OFFSET) * On exit: X(1:OFFSET+LACPTD) contains an approximate * smallest left singular vector of * A(1:OFFSET+LACPTD,1:OFFSET+LACPTD). * * SVLUES (input/output) REAL array, dimension(4) * estimates of singular values. * On entry: SVLUES(1) = sigma_max(A(1:M,1:N)) * SVLUES(2) = sigma_min(A(1:OFFSET,1:OFFSET)) * On exit: SVLUES(1) = sigma_max(A(1:M,1:N)) * SVLUES(2) = sigma_r(B) * SVLUES(3) = sigma_(min(r+1,min(m,n)))(B) * SVLUES(4) = sigma_min(A) * where r = OFFSET+LACPTD and B = A(1:r,1:r) * * MXNM (input/output) FLOATING_DECLARE * On entry: norm of largest column in A(1:OFFSET,1:OFFSET) * On exit: norm of largest column in * A(1:J,1:J) where J = OFFSET+LACPTD * * WORK (workspace) FLOATING_DECLARE array, dimension (LWORK) * * LWORK (input) INTEGER * MAX( 1, N*NB ) if JOB=1, or * MAX( 1, MAX( N, K )*NB ) otherwise. * where NB is the maximum of blocksize used within xGEQRF and * blocksize used within xUNMQR. * * RWORK (workspace) REAL array, dimension ( 2*N ). * * Further Details * =============== * * The matrix Q is represented as a product of elementary reflectors * * Q = H(1) H(2) . . . H(n) * * Each H(i) has the form * * H = I - tau * v * v' * * where tau is a complex scalar, and v is a complex vector with * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i). * * The matrix P is represented in jpvt as follows: If * jpvt(j) = i * then the jth column of P is the ith canonical unit vector. * * ===================================================================== * * .. Parameters .. REAL ZERO, ONE COMPLEX CONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, $ CONE = ( 1.0E+0, 0.0E+0 ) ) * * Indices into the 'svlues' array. * INTEGER IMAX, IBEFOR, IAFTER, IMIN PARAMETER ( IMAX = 1, IBEFOR = 2, IAFTER = 3, IMIN = 4 ) * .. * .. Local Scalars .. INTEGER I, J, PVT, MN, ITEMP, INFO, LASTI REAL TEMP, TEMP2, SMIN, SMINPR, SMAX, SMAXPR COMPLEX AII, SINE, COSINE * .. * .. External Subroutines .. EXTERNAL CLARFG, CLARF, CSWAP, CSCAL, $ CLAIC1, CGEQRF, CUNMQR * .. * .. External Functions .. INTEGER ISAMAX REAL SCNRM2, SLASMX LOGICAL CLAUC1 EXTERNAL ISAMAX, SCNRM2, SLASMX, CLAUC1 * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, CONJG * .. * .. Executable Statements .. * MN = MIN( M, N ) LACPTD = 0 IF( OFFSET.GT.0 ) THEN SMAX = SVLUES( IMAX ) SMIN = SVLUES( IBEFOR ) END IF * * Initialize partial column norms. The first n entries of * work store the exact column norms. * DO 10 I = OFFSET+1,N RWORK( I ) = SCNRM2( M-OFFSET, A( OFFSET+1, I ), 1 ) RWORK( N+I ) = RWORK( I ) 10 CONTINUE * * Compute factorization. * LASTI = MIN( MN, OFFSET+DSRD ) DO 20 I = OFFSET+1, LASTI * * Determine ith pivot column and swap if necessary. * PVT = ( I-1 )+ISAMAX( N-I+1, RWORK( I ), 1 ) IF( PVT.NE.I ) THEN CALL CSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 ) ITEMP = JPVT( PVT ) JPVT( PVT ) = JPVT( I ) JPVT( I ) = ITEMP RWORK( PVT ) = RWORK( I ) RWORK( N+PVT ) = RWORK( N+I ) END IF * * Generate elementary reflector H(i). * IF( I.LT.M ) THEN CALL CLARFG( M-I+1, A( I, I ), A( I+1, I ), 1, $ TAU( I ) ) ELSE CALL CLARFG( 1, A( M, M ), A( M, M ), 1, TAU( M ) ) END IF * * Apply elementary reflector H(I) to the corresponding blocks * of matrices A and C. * AII = A( I, I ) A( I, I ) = CONE IF( I.LT.N ) THEN * * Apply H(I) to A(I:M,I+1:N) from the left. * CALL CLARF( 'Left', M-I+1, N-I, A( I, I ), 1, $ CONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK ) END IF IF( ( JOB.EQ.2 ).AND.( K.GT.0 ) ) THEN * * Apply H(I) to C(I:M,1:K) from the left. * CALL CLARF( 'Left', M-I+1, K, A( I, I ), 1, $ CONJG( TAU( I ) ), C( I, 1 ), LDC, WORK ) ELSE IF( ( JOB.EQ.3 ).AND.( K.GT.0 ) ) THEN * * Apply the transpose of H(I) to C(1:K,I:M) from the right. * CALL CLARF( 'Right', K, M-I+1, A( I, I ), 1, $ TAU( I ), C( 1, I ), LDC, WORK ) END IF A( I, I ) = AII * * Update partial column norms. * IF( I.LT.LASTI ) THEN DO 30 J = I+1, N IF( RWORK( J ).NE.ZERO ) THEN TEMP = ONE-( ABS( A( I, J ) )/RWORK( J ) )**2 TEMP = MAX( TEMP, ZERO ) TEMP2 = ONE+0.05*TEMP*( RWORK( J )/RWORK( N+J ) )**2 IF( TEMP2.EQ.ONE ) THEN RWORK( J ) = SCNRM2( M-I, A( I+1, J ), 1 ) RWORK( N+J ) = RWORK( J ) ELSE RWORK( J ) = RWORK( J )*SQRT( TEMP ) END IF END IF 30 CONTINUE END IF * * Check new column for independence. * IF( I.EQ.1 ) THEN MXNM = ABS( A( 1, 1 ) ) SMIN = MXNM SMAX = MXNM X( 1 ) = CONE IF( MXNM.GT.ZERO ) THEN LACPTD = 1 ELSE SVLUES( IAFTER ) = SMIN GOTO 50 END IF ELSE SMAXPR = SLASMX( I )*MXNM IF( CLAUC1( I, X, SMIN, A( 1, I ), A( I, I ), $ SMAXPR*IRCOND ) ) THEN * * Column accepted. * SMAX = SMAXPR LACPTD = LACPTD + 1 ELSE * * Column rejected. * GOTO 50 END IF END IF 20 CONTINUE 50 SVLUES( IMAX ) = SMAX SVLUES( IBEFOR ) = SMIN IF( LACPTD.EQ.DSRD ) THEN * * DSRD independent columns have been found. * SVLUES( IAFTER ) = SMIN SVLUES( IMIN ) = SMIN ELSE * * All remaining columns rejected. * I = OFFSET + LACPTD + 1 IF( I.LT.MN ) THEN * * Factor remaining columns. * CALL CGEQRF( M-I, N-I, A( I+1, I+1 ), LDA, TAU( I+1 ), $ WORK, LWORK, INFO ) * * Apply the transformations computed in CGEQRF to the * corresponding part of matrix C. * IF( ( JOB.EQ.2 ).AND.( K.GT.0 ) ) THEN * * Apply them to C(I+1:M,1:K) from the left. * CALL CUNMQR( 'Left', 'Conjugate Transpose', $ M-I, K, MN-I, $ A( I+1, I+1 ), LDA, TAU( I+1 ), $ C( I+1, 1 ), LDC, WORK, LWORK, INFO ) ELSE IF( ( JOB.EQ.3 ).AND.( K.GT.0 ) ) THEN * * Apply the transpose of them to C(1:K,I+1:M) from the * right. * CALL CUNMQR( 'Right', 'No Transpose', K, M-I, MN-I, $ A( I+1, I+1 ), LDA, TAU( I+1 ), $ C( 1, I+1 ), LDC, WORK, LWORK, INFO ) END IF END IF * * Use incremental condition estimation to get an estimate * of the smallest singular value. * DO 60 I = MAX( 2, OFFSET+LACPTD+1 ), MN CALL CLAIC1( 2, I-1, X, SMIN, A( 1, I ), A( I, I ), $ SMINPR, SINE, COSINE ) CALL CSCAL( I-1, SINE, X, 1 ) X( I ) = COSINE SMIN = SMINPR IF( I.EQ.OFFSET+LACPTD+1 ) THEN SVLUES( IAFTER ) = SMIN END IF 60 CONTINUE SVLUES( IMIN ) = SMIN END IF RETURN * * End of CGEQPC * END SHAR_EOF fi # end of overwriting check if test -f 'cgeqpw.f' then echo shar: will not over-write existing file "'cgeqpw.f'" else cat << SHAR_EOF > 'cgeqpw.f' SUBROUTINE CGEQPW( M, LWSIZE, NB, OFFSET, LACPTD, A, LDA, JPVT, $ IRCOND, X, SMIN, MXNM, TAU, WORK, RWORK ) * * This code is part of a release of the package for computing * rank-revealing QR Factorizations written by: * ================================================================== * Christian H. Bischof and Gregorio Quintana-Orti * Math. and Comp. Sci. Div. Departamento de Informatica * Argonne National Lab. Universidad Jaime I * Argonne, IL 60439 Campus P. Roja, 12071 Castellon * USA Spain * bischof@mcs.anl.gov gquintan@inf.uji.es * ================================================================== * $Revision: 1.42 $ * $Date: 96/12/30 16:59:38 $ * * .. Scalar Arguments .. INTEGER M, LWSIZE, NB, OFFSET, LACPTD, LDA REAL IRCOND, SMIN, MXNM * .. * .. Array Arguments .. INTEGER JPVT( * ) COMPLEX A( LDA, * ), TAU( * ), X( * ), WORK( * ) REAL RWORK( * ) * * * Purpose * ======= * * CGEQPW applies one block step of the Householder QR * factorization algorithm with restricted pivoting. It is called * by CGEQPB to factorize a window of the matrix. * * Let A be the partial QR factorization of an M by (OFFSET+LWSIZE) * matrix C, i.e. we have computed an orthogonal matrix Q1 and a * permutation matrix P1 such that * C * P1 = Q1 * A * and A(:,1:OFFSET) is upper triangular. Let us denote A(:,1:OFFSET) * by B. Then in addition let * X be an approximate smallest left singular vector of B in the sense * that * sigma_min(B) ~ twonorm(B'*X) = SMIN * and * sigma_max(B) ~ ((offset)**(1./3.))*MXNM = SMAX * with * cond_no(B) ~ SMAX/SMIN <= 1/IRCOND * * Then CGEQP2 tries to identify NB columns in * A(:,OFFSET+1:OFFSET+LWSIZE) such that * cond_no([B,D]) < 1/IRCOND * where D are the KB columns of A(:,OFFSET+1:OFFSET+LWSIZE) that were * considered independent with respect to the threshold 1/IRCOND. * * On exit, * C * P2 = Q2 * A * is again a partial QR factorization of C, but columns * OFFSET+1:OFFSET+LACPTD of A have been reduced via * a series of elementary reflectors to upper * trapezoidal form. Further * sigma_min(A(:,1:OFFSET+LACPTD)) * ~ twonorm(A(:,1:OFFSET+LACPTD)'*x) = SMIN * and * sigma_max(A(:,1:OFFSET+LACPTD)) ~ sqrt(OFFSET+LACPTD)*MXNM = SMAX * with * cond_no(A(:,1:OFFSET+LACPTD)) * ~ SMAX/SMIN <= 1/IRCOND. * * In the ideal case, LACPTD = NB, that is, * we found NB independent columns in the window consisting of * the first LWSIZE columns of A. * * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. * * LWSIZE (input) INTEGER * The size of the pivot window in A. * * NB (input) INTEGER * The number of independent columns one would like to identify. * This equals the desired blocksize in CGEQPB. * * OFFSET (input) INTEGER * The number of rows and columns of A that need not be updated. * * LACPTD (output) INTEGER * The number of columns in A(:,OFFSET+LWSIZE) that were * accepted as linearly independent. * * A (input/output) COMPLEX array, dimension (LDA,OFFSET+LWSIZE) * On entry, the upper triangle of A(:,1:OFFSET) contains the * partially completed triangular factor R; the elements below * the diagonal, with the array TAU, represent the matrix Q1 as * a product of elementary reflectors. * On exit, the upper triangle of A(:,OFFSET+LACPTD) contains * the partially completed upper triangular factor R; the * elements below the diagonal, with the array TAU, represent * the matrix Q2 as a product of elementary reflectors. * A(OFFSET:M,LACPTD+1:LWSIZE) has been updated by the product * of these elementary reflectors. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= M. * * JPVT (input/output) INTEGER array, dimension (OFFSET+LWSIZE) * On entry and exit, jpvt(i) = k if the i-th column * of A was the k-th column of C. * * IRCOND (input) REAL * 1/IRCOND is the threshold for the condition number. * * X (input/output) COMPLEX array, dimension (OFFSET+NB) * On entry, X(1:OFFSET) is an approximate left nullvector of * the upper triangle of A(1:OFFSET,1:OFFSET). * On exit, X(1:OFFSET+LACPTD) is an approximate left * nullvector of the matrix in the upper triangle of * A(1:OFFSET+LACPTD,1:OFFSET+LACPTD). * * SMIN (input/output) REAL * On entry, SMIN is an estimate for the smallest singular * value of the upper triangle of A(1:OFFSET,1:OFFSET). * On exit, SMIN is an estimate for the smallest singular * value of the matrix in the upper triangle of * A(1:OFFSET+LACPTD,1:OFFSET+LACPTD). * * MXNM (input) FLOATING_DECLARE * The norm of the largest column in matrix A. * * TAU (output) COMPLEX array, dimension (OFFSET+LWSIZE) * On exit, TAU(1:OFFSET+LACPTD) contains details of * the matrix Q2. * * WORK (workspace) COMPLEX array, dimension (LWSIZE) * * RWORK (workspace) REAL array, dimension (2*LWSIZE) * * ================================================================ * * .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. INTEGER I, K, I1, PVTIDX, LASTK REAL TEMP, TEMP2, SMAX COMPLEX GAMMA, AKK * .. * .. External Subroutines .. EXTERNAL SCNRM2, CSCAL, CSWAP, CLARFG, $ CLARF, ISAMAX, CLAUC1, SLASMX INTEGER ISAMAX REAL SCNRM2, SLASMX LOGICAL CLAUC1 * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT, REAL, CMPLX, CONJG * .. * .. Executable Statements .. * * Initialize partial column norms (stored in the first LWSIZE * entries of WORK) and exact column norms (stored in the second * LWSIZE entries of WORK) for the first batch of columns. * DO 10 I = 1,LWSIZE RWORK( I ) = SCNRM2( M-OFFSET, A( OFFSET+1, OFFSET+I ), 1 ) RWORK( LWSIZE+I ) = RWORK( I ) 10 CONTINUE * * ************* * * Main loop * * ************* * LASTK = MIN( M, OFFSET+LWSIZE ) LACPTD = 0 1000 IF( LACPTD.EQ.NB ) GOTO 2000 * * Determine pivot candidate. * ========================= PVTIDX = OFFSET + LACPTD + $ ISAMAX( LWSIZE-LACPTD, RWORK( LACPTD+1 ), 1 ) K = OFFSET + LACPTD + 1 * * Exchange current column and pivot column. * IF( PVTIDX.NE.K ) THEN CALL CSWAP( M, A( 1, PVTIDX ), 1, A( 1, K ), 1 ) I1 = JPVT( PVTIDX ) JPVT( PVTIDX ) = JPVT( K ) JPVT( K ) = I1 TEMP = RWORK( PVTIDX-OFFSET ) RWORK( PVTIDX-OFFSET ) = RWORK( K-OFFSET ) RWORK( K-OFFSET ) = TEMP TEMP = RWORK( PVTIDX-OFFSET+LWSIZE ) RWORK( PVTIDX-OFFSET+LWSIZE ) = RWORK( K+LWSIZE-OFFSET ) RWORK( K+LWSIZE-OFFSET ) = TEMP END IF * * Determine (offset+lacptd+1)st diagonal element * GAMMA of matrix A should elementary reflector be applied. * TEMP = REAL( A( K, K ) ) IF( TEMP.EQ.ZERO ) THEN GAMMA = -RWORK( K-OFFSET ) ELSE GAMMA = -SIGN( RWORK( K-OFFSET ), TEMP ) END IF * * Update estimate for largest singular value. * SMAX = SLASMX( K )*MXNM * * Is candidate pivot column acceptable ? * ===================================== IF( CLAUC1( K, X, SMIN, A( 1, K ), GAMMA, SMAX*IRCOND ) ) $ THEN * * Pivot candidate was accepted. * ============================ * LACPTD = LACPTD + 1 * * Generate Householder vector. * IF( K.LT.M ) THEN CALL CLARFG( M-K+1, A( K, K ), A( K+1, K ), 1, $ TAU( K ) ) ELSE CALL CLARFG( 1, A( M, K), A( M, K ), 1, TAU ( K ) ) END IF * * Apply Householder reflection to A(k:m,k+1:lwsize). * IF( LACPTD.LT.LWSIZE ) THEN AKK = A( K, K ) A( K, K ) = CMPLX( ONE ) CALL CLARF( 'Left', M-K+1, LWSIZE-LACPTD, $ A( K, K ), 1, CONJG( TAU( K ) ), $ A( K, K+1 ), LDA, WORK ) A( K, K ) = AKK END IF * * Update partial column norms. * IF( K.LT.LASTK ) THEN DO 20 I = LACPTD+1,LWSIZE IF( RWORK( I ).NE.ZERO ) THEN TEMP = ONE- $ ( ABS( A( K, OFFSET+I ) )/RWORK( I ) )**2 TEMP = MAX( TEMP, ZERO ) TEMP2 = ONE+ 0.05*TEMP* $ ( RWORK( I )/RWORK( I+LWSIZE ) )**2 IF( TEMP2.EQ.ONE ) THEN RWORK( I ) = SCNRM2( M-K, $ A( K+1, OFFSET+I ), 1 ) RWORK( I+LWSIZE ) = RWORK( I ) ELSE RWORK( I ) = RWORK( I )*SQRT( TEMP ) END IF END IF 20 CONTINUE END IF ELSE * * Reject all remaining columns in pivot window. * ============================================ * GOTO 2000 END IF * * End while. * GOTO 1000 2000 CONTINUE RETURN * * End of CGEQPW * END SHAR_EOF fi # end of overwriting check if test -f 'cgeqpx.f' then echo shar: will not over-write existing file "'cgeqpx.f'" else cat << SHAR_EOF > 'cgeqpx.f' SUBROUTINE CGEQPX( JOB, M, N, K, A, LDA, C, LDC, JPVT, IRCOND, $ ORCOND, RANK, SVLUES, WORK, LWORK, RWORK, $ INFO ) * * This code is part of a release of the package for computing * rank-revealing QR Factorizations written by: * ================================================================== * Christian H. Bischof and Gregorio Quintana-Orti * Math. and Comp. Sci. Div. Departamento de Informatica * Argonne National Lab. Universidad Jaime I * Argonne, IL 60439 Campus P. Roja, 12071 Castellon * USA Spain * bischof@mcs.anl.gov gquintan@inf.uji.es * ================================================================== * $Revision: 1.42 $ * $Date: 96/12/30 16:59:38 $ * * .. Scalar Arguments .. INTEGER JOB, M, N, K, LDA, LDC, RANK, LWORK, INFO REAL IRCOND, ORCOND * .. * .. Array Arguments .. INTEGER JPVT( * ) COMPLEX A( LDA, * ), C( LDC, * ), WORK( * ) REAL SVLUES( 4 ), RWORK( * ) * .. * * Purpose * ======= * * CGEQPX computes a QR factorization * A*P = Q*[ R11 R12 ] * [ 0 R22 ] * of a real m by n matrix A. The permutation P is * chosen with the goal to reveal the rank of A by a * suitably dimensioned trailing submatrix R22 with norm(R22) * being small. * * Based on methods related to Chandrasekaran&Ipsen's algorithms. * * Arguments * ========= * * JOB (input) INTEGER * The job to do: * = 1: The transformations needed in the * triangularization are only applied to matrix A. * Thus, matrix C is not updated. * = 2: The same transformations needed in the * triangularization of matrix A are applied to * matrix C from the left. * That is, if Q'*A*P=R, then C := Q'*C. * In this case, matrix C is m-by-k. * = 3: The transpose of the transformations needed * in the triangularization of matrix A are applied * to matrix C from the right. * That is, if Q'*A*P=R, then C := C*Q. * In this case, matrix C is k-by-m. * In these three cases, the permutations are always saved * into vector JPVT. * * M (input) INTEGER * The number of rows of matrices A. M >= 0. * If JOB=2, M is the number of rows of matrix C. * If JOB=3, M is the number of columns of matrix C. * * N (input) INTEGER * The number of columns of matrix A. N >= 0. * * K (input) INTEGER * It defines the dimension of matrix C. K >= 0. * If JOB=2, K is the number of columns of matrix C. * If JOB=3, K is the number of rows of matrix C. * * A (input/output) COMPLEX array, dimension (LDA,N) * On entry, the m by n matrix A. * On exit, the upper triangle of the array contains the * min(m,n) by n upper trapezoidal matrix R; the lower triangle * array is filled with zeros. * * LDA (input) INTEGER * The leading dimension of array A. LDA >= MAX(1,M). * * C (input/output) COMPLEX array, dimension * ( LDC, K ) if JOB=2. * ( LDC, M ) if JOB=3. * If argument JOB asks, all the transformations * applied to matrix A are also applied to matrix C. * * LDC (input) INTEGER * The leading dimension of array C. * If JOB=2, then LDC >= MAX(1,M). * If JOB=3, then LDC >= MAX(1,K). * * JPVT (output) INTEGER array, dimension (N) * JPVT(I) = J <==> Column J of A has been permuted into * position I in AP. * JPVT(1:RANK) contains the indices of the columns considered * linearly independent. * JPVT(RANK+1:N) contains the indices of the columns considered * linearly dependent from the previous ones. * * IRCOND (input) FLOATING_DECLARE * 1/IRCOND specifies an upper bound on the condition number * of R11. If IRCOND == 0, IRCOND = machine precision is chosen * as default. IRCOND must be >= 0. * * ORCOND (output) FLOATING_DECLARE * 1/ORCOND is an estimate for the condition number of R11. * * RANK (output) INTEGER * RANK is an estimate for the numerical rank of A with respect * to the threshold 1/IRCOND in the sense that * RANK = arg_max(cond_no(R(1:r,1:r))<1/IRCOND) * * SVLUES (output) REAL array, dimension (4) * On exit, SVLUES contains estimates of some of the * singular values of the triangular factor R. * SVLUES(1): largest singular value of R(1:RANK,1:RANK) * SVLUES(2): smallest singular value of R(1:RANK,1:RANK) * SVLUES(3): smallest singular value of R(1:RANK+1,1:RANK+1) * SVLUES(4): smallest singular value of R * If the triangular factorization is a rank-revealing one * (which will be the case if the leading columns were well- * conditioned), then SVLUES(1) will also be an estimate * for the largest singular value of A, SVLUES(2) and SVLUES(3) * will be estimates for the RANK-th and (RANK+1)-st singular * value of A, and SVLUES(4) wil be an estimate for the * smallest singular value of A. * By examining these values, one can confirm that the rank is * well defined with respect to the threshold chosen. * * WORK (workspace) COMPLEX array, dimension (LWORK) * On exit: work(1) is the size of the storage array needed * for optimal performance * * LWORK (input) INTEGER * The dimension of array WORK. * If JOB=1: * The unblocked strategy requires that: * LWORK >= 2*MN+N. * The block algorithm requires that: * LWORK >= 2*MN+N*NB. * If JOB<>1: * The unblocked strategy requires that: * LWORK >= 2*MN+MAX(K,N). * The block algorithm requires that: * LWORK >= 2*MN+NB*NB+NB*MAX(K,N). * Where MN = min(M,N) and NB is the block size for this * environment. * In both cases, the minimum required workspace is the * one for the unblocked strategy. * * RWORK (workspace) REAL array, dimension ( 2*N ) * * INFO (output) INTEGER * = 0: Successful exit. * < 0: If INFO = -i, the i-th argument had an illegal value * > 0: Problems in the computation of the rank. * 1: Exceeded the allowed maximum number of steps. * 2: Rank not well defined. * In adition, vector SVLUES tell if rank is not well defined. * * ===================================================================== * * .. Parameters .. REAL ZERO PARAMETER ( ZERO = 0.0E+0 ) * .. * .. External Subroutines .. EXTERNAL XERBLA, CGEQPB, CTRQPX * .. * .. Intrinsic Functions .. INTRINSIC MIN, MAX, REAL * .. * .. Local Scalars .. REAL WSIZE INTEGER MN, WKMIN * .. * .. Executable Statements .. * MN = MIN( M, N ) IF( JOB.EQ.1 ) THEN WKMIN = 2*MN+N ELSE WKMIN = 2*MN+MAX(K,N) END IF * * Test input arguments * ==================== * INFO = 0 IF( ( JOB.LT.1 ).OR.( JOB.GT.3 ) ) THEN INFO = -1 ELSE IF( M.LT.0 ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( K.LT.0 ) THEN INFO = -4 ELSE IF( LDA.LT.MAX(1,M) ) THEN INFO = -6 ELSE IF( ( ( JOB.EQ.1 ).AND.( LDC.LT.1 ) ).OR. $ ( ( JOB.EQ.2 ).AND.( LDC.LT.MAX( 1, M ) ) ).OR. $ ( ( JOB.EQ.3 ).AND.( LDC.LT.MAX( 1, K ) ) ) ) THEN INFO = -8 ELSE IF( IRCOND.LT.ZERO ) THEN INFO = -10 ELSE IF( LWORK.LT.MAX( 1, WKMIN ) ) THEN INFO = -15 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CGEQPX',-INFO ) RETURN END IF * * Preprocessing * ============= * CALL CGEQPB( JOB, M, N, K, A, LDA, C, LDC, JPVT, IRCOND, $ ORCOND, RANK, SVLUES, WORK, LWORK, RWORK, INFO ) WSIZE = REAL( WORK( 1 ) ) * * Postprocessing * ============== * IF( RANK.GT.0 ) THEN CALL CTRQPX( JOB, M, N, K, A, LDA, C, LDC, JPVT, IRCOND, $ ORCOND, RANK, SVLUES, WORK, RWORK, INFO ) END IF * WORK( 1 ) = WSIZE RETURN * * End of CGEQPX * END SHAR_EOF fi # end of overwriting check if test -f 'cgeqpy.f' then echo shar: will not over-write existing file "'cgeqpy.f'" else cat << SHAR_EOF > 'cgeqpy.f' SUBROUTINE CGEQPY( JOB, M, N, K, A, LDA, C, LDC, JPVT, IRCOND, $ ORCOND, RANK, SVLUES, WORK, LWORK, RWORK, $ INFO ) * * This code is part of a release of the package for computing * rank-revealing QR Factorizations written by: * ================================================================== * Christian H. Bischof and Gregorio Quintana-Orti * Math. and Comp. Sci. Div. Departamento de Informatica * Argonne National Lab. Universidad Jaime I * Argonne, IL 60439 Campus P. Roja, 12071 Castellon * USA Spain * bischof@mcs.anl.gov gquintan@inf.uji.es * ================================================================== * $Revision: 1.42 $ * $Date: 96/12/30 16:59:40 $ * * .. Scalar Arguments .. INTEGER JOB, M, N, K, LDA, LDC, RANK, LWORK, INFO REAL IRCOND, ORCOND * .. * .. Array Arguments .. INTEGER JPVT( * ) COMPLEX A( LDA, * ), C( LDC, * ), WORK( * ) REAL SVLUES( 4 ), RWORK( * ) * .. * * Purpose * ======= * * CGEQPY computes a QR factorization * A*P = Q*[ R11 R12 ] * [ 0 R22 ] * of a real m by n matrix A. The permutation P is * chosen with the goal to reveal the rank of A by a * suitably dimensioned trailing submatrix R22 with norm(R22) * being small. * * Based on Pan&Tang's algorithm number 3. * * Arguments * ========= * * JOB (input) INTEGER * The job to do: * = 1: The transformations needed in the * triangularization are only applied to matrix A. * Thus, matrix C is not updated. * = 2: The same transformations needed in the * triangularization of matrix A are applied to * matrix C from the left. * That is, if Q'*A*P=R, then C := Q'*C. * In this case, matrix C is m-by-k. * = 3: The transpose of the transformations needed * in the triangularization of matrix A are applied * to matrix C from the right. * That is, if Q'*A*P=R, then C := C*Q. * In this case, matrix C is k-by-m. * In these three cases, the permutations are always saved * into vector JPVT. * * M (input) INTEGER * The number of rows of matrices A. M >= 0. * If JOB=2, M is the number of rows of matrix C. * If JOB=3, M is the number of columns of matrix C. * * N (input) INTEGER * The number of columns of matrix A. N >= 0. * * K (input) INTEGER * It defines the dimension of matrix C. K >= 0. * If JOB=2, K is the number of columns of matrix C. * If JOB=3, K is the number of rows of matrix C. * * A (input/output) COMPLEX array, dimension (LDA,N) * On entry, the m by n matrix A. * On exit, the upper triangle of the array contains the * min(m,n) by n upper trapezoidal matrix R; the lower triangle * array is filled with zeros. * * LDA (input) INTEGER * The leading dimension of array A. LDA >= MAX(1,M). * * C (input/output) COMPLEX array, dimension * ( LDC, K ) if JOB=2. * ( LDC, M ) if JOB=3. * If argument JOB asks, all the transformations * applied to matrix A are also applied to matrix C. * * LDC (input) INTEGER * The leading dimension of array C. * If JOB=2, then LDC >= MAX(1,M). * If JOB=3, then LDC >= MAX(1,K). * * JPVT (output) INTEGER array, dimension (N) * JPVT(I) = J <==> Column J of A has been permuted into * position I in AP. * JPVT(1:RANK) contains the indices of the columns considered * linearly independent. * JPVT(RANK+1:N) contains the indices of the columns considered * linearly dependent from the previous ones. * * IRCOND (input) FLOATING_DECLARE * 1/IRCOND specifies an upper bound on the condition number * of R11. If IRCOND == 0, IRCOND = machine precision is chosen * as default. IRCOND must be >= 0. * * ORCOND (output) FLOATING_DECLARE * 1/ORCOND is an estimate for the condition number of R11. * * RANK (output) INTEGER * RANK is an estimate for the numerical rank of A with respect * to the threshold 1/IRCOND in the sense that * RANK = arg_max(cond_no(R(1:r,1:r))<1/IRCOND) * * SVLUES (output) REAL array, dimension (4) * On exit, SVLUES contains estimates of some of the * singular values of the triangular factor R. * SVLUES(1): largest singular value of R(1:RANK,1:RANK) * SVLUES(2): smallest singular value of R(1:RANK,1:RANK) * SVLUES(3): smallest singular value of R(1:RANK+1,1:RANK+1) * SVLUES(4): smallest singular value of R * If the triangular factorization is a rank-revealing one * (which will be the case if the leading columns were well- * conditioned), then SVLUES(1) will also be an estimate * for the largest singular value of A, SVLUES(2) and SVLUES(3) * will be estimates for the RANK-th and (RANK+1)-st singular * value of A, and SVLUES(4) wil be an estimate for the * smallest singular value of A. * By examining these values, one can confirm that the rank is * well defined with respect to the threshold chosen. * * WORK (workspace) COMPLEX array, dimension (LWORK) * On exit: work(1) is the size of the storage array needed * for optimal performance * * LWORK (input) INTEGER * The dimension of array WORK. * If JOB=1: * The unblocked strategy requires that: * LWORK >= 2*MN+N. * The block algorithm requires that: * LWORK >= 2*MN+N*NB. * If JOB<>1: * The unblocked strategy requires that: * LWORK >= 2*MN+MAX(K,N). * The block algorithm requires that: * LWORK >= 2*MN+NB*NB+NB*MAX(K,N). * Where MN = min(M,N) and NB is the block size for this * environment. * In both cases, the minimum required workspace is the * one for the unblocked strategy. * * RWORK (workspace) REAL array, dimension ( 2*N ) * * INFO (output) INTEGER * = 0: Successful exit. * < 0: If INFO = -i, the i-th argument had an illegal value * > 0: Problems in the computation of the rank. * 1: Exceeded the allowed maximum number of steps. * 2: Rank not well defined. * In adition, vector SVLUES tell if rank is not well defined. * * ===================================================================== * * .. Parameters .. REAL ZERO PARAMETER ( ZERO = 0.0E+0 ) * .. * .. External Subroutines .. EXTERNAL XERBLA, CGEQPB, CTRQPY * .. * .. Intrinsic Functions .. INTRINSIC MIN, MAX, REAL * .. * .. Local Scalars .. REAL WSIZE INTEGER MN, WKMIN * .. * .. Executable Statements .. * MN = MIN( M, N ) IF( JOB.EQ.1 ) THEN WKMIN = 2*MN+N ELSE WKMIN = 2*MN+MAX(K,N) END IF * * Test input arguments * ==================== * INFO = 0 IF( ( JOB.LT.1 ).OR.( JOB.GT.3 ) ) THEN INFO = -1 ELSE IF( M.LT.0 ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( K.LT.0 ) THEN INFO = -4 ELSE IF( LDA.LT.MAX(1,M) ) THEN INFO = -6 ELSE IF( ( ( JOB.EQ.1 ).AND.( LDC.LT.1 ) ).OR. $ ( ( JOB.EQ.2 ).AND.( LDC.LT.MAX( 1, M ) ) ).OR. $ ( ( JOB.EQ.3 ).AND.( LDC.LT.MAX( 1, K ) ) ) ) THEN INFO = -8 ELSE IF( IRCOND.LT.ZERO ) THEN INFO = -10 ELSE IF( LWORK.LT.MAX( 1, WKMIN ) ) THEN INFO = -15 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CGEQPY',-INFO ) RETURN END IF * * Preprocessing * ============= * CALL CGEQPB( JOB, M, N, K, A, LDA, C, LDC, JPVT, IRCOND, $ ORCOND, RANK, SVLUES, WORK, LWORK, RWORK, INFO ) WSIZE = REAL( WORK( 1 ) ) * * Postprocessing * ============== * IF( RANK.GT.0 ) THEN CALL CTRQPY( JOB, M, N, K, A, LDA, C, LDC, JPVT, IRCOND, $ ORCOND, RANK, SVLUES, WORK, RWORK, INFO ) END IF * WORK( 1 ) = WSIZE RETURN * * End of CGEQPY * END SHAR_EOF fi # end of overwriting check if test -f 'clasmx.f' then echo shar: will not over-write existing file "'clasmx.f'" else cat << SHAR_EOF > 'clasmx.f' REAL FUNCTION SLASMX( I ) INTEGER I * REAL OTHIRD PARAMETER ( OTHIRD = 1.0E+0/3.0E+0 ) INTRINSIC REAL SLASMX = REAL( I )**OTHIRD RETURN END SHAR_EOF fi # end of overwriting check if test -f 'clauc1.f' then echo shar: will not over-write existing file "'clauc1.f'" else cat << SHAR_EOF > 'clauc1.f' LOGICAL FUNCTION CLAUC1( K, X, SMIN, W, GAMMA, THRESH ) * * This code is part of a release of the package for computing * rank-revealing QR Factorizations written by: * ================================================================== * Christian H. Bischof and Gregorio Quintana-Orti * Math. and Comp. Sci. Div. Departamento de Informatica * Argonne National Lab. Universidad Jaime I * Argonne, IL 60439 Campus P. Roja, 12071 Castellon * USA Spain * bischof@mcs.anl.gov gquintan@inf.uji.es * ================================================================== * $Revision: 1.42 $ * $Date: 96/12/30 16:59:41 $ * * .. Scalar Arguments .. INTEGER K REAL SMIN, THRESH COMPLEX GAMMA * .. * .. Array Arguments .. COMPLEX W( * ), X( * ) * .. * * Purpose * ======= * * PREC_LAUC1 applies incremental condition estimation to determine whether * the K-th column of A, stored in vector W, would be acceptable as a pivot * column with respect to the threshold THRESH. * * Arguments * ========= * * K (input) INTEGER * Length of vector X. * * X (input/output) COMPLEX array, dimension ( K ) * On entry, X(1:K-1) contains an approximate smallest left singular * vector of the upper triangle of A(1:k-1,1:k-1). * On exit, if CLAUC1 == .TRUE., X contains an approximate * smallest left singular vector of the upper triangle of A(1:k,1:k); * if CLAUC1 == .FALSE., X is unchanged. * * SMIN (input/output) REAL * On entry, an estimate for the smallest singular value of the * upper triangle of A(1:k-1,1:k-1). * On exit, if CLAUC1 == .TRUE., SMIN is an estimate of the * smallest singular value of the upper triangle of A(1:k,1:k); * if CLAUC1 == .FALSE., SMIN is unchanged. * * W (input) FLOATING_DECLARE array, dimension ( K-1 ) * The K-th column of matrix A excluding the diagonal element. * * GAMMA (input) COMPLEX * Diagonal entry in k-th column of A if column k were to * be accepted. * * THRESH (input) REAL * If the approximate smallest singular value for A(1:K,1:K) * is smaller than THRESH, the kth column is rejected. * * (CLAUC1) (output) LOGICAL * If the k-th column of A is found acceptable, CLAUC1 * returns .TRUE., otherwise it returns .FALSE. * * ===================================================================== * * .. Local Scalars .. REAL SMINPR COMPLEX SINE, COSINE * .. * .. External Subroutines .. EXTERNAL CLAIC1, CSCAL * .. * .. Intrinsic Functions .. INTRINSIC ABS * .. * .. Executable Statements .. * * * Try to use diagonal element as condition estimator * IF( THRESH.GT.ABS( GAMMA ) ) THEN CLAUC1 = .FALSE. RETURN END IF * * Use incremental condition estimation to determine an estimate * SMINPR and an approximate singular vector [SINE*X,COSINE]' * for A(K,K). * CALL CLAIC1( 2, K-1, X, SMIN, W, GAMMA, SMINPR, $ SINE, COSINE ) IF( THRESH.GT.SMINPR ) THEN CLAUC1 = .FALSE. ELSE CALL CSCAL( K-1, SINE, X, 1 ) X( K ) = COSINE SMIN = SMINPR CLAUC1 = .TRUE. END IF RETURN * * End of CLAUC1 * END SHAR_EOF fi # end of overwriting check if test -f 'cmylap.f' then echo shar: will not over-write existing file "'cmylap.f'" else cat << SHAR_EOF > 'cmylap.f' ********************************************************************* SUBROUTINE CGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO ) * * -- LAPACK auxiliary routine (version 1.0b) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER JPVT( * ) REAL RWORK( * ) COMPLEX A( LDA, * ), TAU( * ), WORK( * ) * .. * * Purpose * ======= * * CGEQPF computes a QR factorization with column pivoting of a * complex m by n matrix A: A*P = Q*R * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0 * * A (input/output) COMPLEX array, dimension (LDA,N) * On entry, the m by n matrix A. * On exit, the upper triangle of the array contains the * min(m,n) by n upper triangular matrix R; the elements * below the diagonal, together with the array TAU, * represent the orthogonal matrix Q as a product of * min(m,n) elementary reflectors. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * JPVT (input/output) INTEGER array, dimension (N) * on entry: If JPVT(I) <> 0, column I of A is permuted * to the front of AP (a leading column) * IF JPVT(I) == 0, column I of A is a free column. * on exit: If JPVT(I) = K, then the Ith column of AP * was the Kth column of A. * * TAU (output) COMPLEX array, dimension (min(M,N)) * Stores further details of * the orthogonal matrix Q (see A). * * WORK (workspace) COMPLEX array, dimension (N) * * RWORK (workspace) REAL array, dimension (2*N) * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * Further Details * =============== * * The matrix Q is represented as a product of elementary reflectors * * Q = H(1) H(2) . . . H(n) * * Each H(i) has the form * * H = I - tau * v * v' * * where tau is a complex scalar, and v is a complex vector with * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i). * * The matrix P is represented in jpvt as follows: If * jpvt(j) = i * then the jth column of P is the ith canonical unit vector. * * ===================================================================== * * .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. INTEGER I, ITEMP, J, MA, MN, PVT REAL TEMP, TEMP2 COMPLEX AII * .. * .. External Subroutines .. EXTERNAL CGEQR2, CLARF, CLARFG, CSWAP, CUNM2R, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, CMPLX, CONJG, MAX, MIN, SQRT * .. * .. External Functions .. INTEGER ISAMAX REAL SCNRM2 EXTERNAL ISAMAX, SCNRM2 * .. * .. Executable Statements .. * * Test the input arguments * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CGEQPF', -INFO ) RETURN END IF * * Quick return if possible * IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * MN = MIN( M, N ) * * Move initial columns up front * ITEMP = 1 DO 10 I = 1, N IF( JPVT( I ).NE.0 ) THEN IF( I.NE.ITEMP ) THEN CALL CSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 ) JPVT( I ) = JPVT( ITEMP ) JPVT( ITEMP ) = I ELSE JPVT( I ) = I END IF ITEMP = ITEMP + 1 ELSE JPVT( I ) = I END IF 10 CONTINUE ITEMP = ITEMP - 1 * * Compute the QR factorization and update remaining columns * IF( ITEMP.GT.0 ) THEN MA = MIN( ITEMP, M ) CALL CGEQR2( M, MA, A, LDA, TAU, WORK, INFO ) IF( MA.LT.N ) THEN CALL CUNM2R( 'Left', 'Conjugate transpose', M, N-MA, MA, A, $ LDA, TAU, A( 1, MA+1 ), LDA, WORK, INFO ) END IF END IF * IF( ITEMP.LT.MN ) THEN * * Initialize partial column norms. The first n entries of * work store the exact column norms. * DO 20 I = ITEMP + 1, N RWORK( I ) = SCNRM2( M-ITEMP, A( ITEMP+1, I ), 1 ) RWORK( N+I ) = RWORK( I ) 20 CONTINUE * * Compute factorization * DO 40 I = ITEMP + 1, MN * * Determine ith pivot column and swap if necessary * PVT = ( I-1 ) + ISAMAX( N-I+1, RWORK( I ), 1 ) * IF( PVT.NE.I ) THEN CALL CSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 ) ITEMP = JPVT( PVT ) JPVT( PVT ) = JPVT( I ) JPVT( I ) = ITEMP RWORK( PVT ) = RWORK( I ) RWORK( N+PVT ) = RWORK( N+I ) END IF * * Generate elementary reflector H(i) * AII = A( I, I ) CALL CLARFG( M-I+1, AII, A( MIN( I+1, M ), I ), 1, $ TAU( I ) ) A( I, I ) = AII * IF( I.LT.N ) THEN * * Apply H(i) to A(i:m,i+1:n) from the left * AII = A( I, I ) A( I, I ) = CMPLX( ONE ) CALL CLARF( 'Left', M-I+1, N-I, A( I, I ), 1, $ CONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK ) A( I, I ) = AII END IF * * Update partial column norms * DO 30 J = I + 1, N IF( RWORK( J ).NE.ZERO ) THEN TEMP = ONE - ( ABS( A( I, J ) ) / RWORK( J ) )**2 TEMP = MAX( TEMP, ZERO ) TEMP2 = ONE + 0.05*TEMP*( RWORK( J ) / RWORK( N+J ) ) $ **2 IF( TEMP2.EQ.ONE ) THEN IF( M-I.GT.0 ) THEN RWORK( J ) = SCNRM2( M-I, A( I+1, J ), 1 ) RWORK( N+J ) = RWORK( J ) ELSE RWORK( J ) = ZERO RWORK( N+J ) = ZERO END IF ELSE RWORK( J ) = RWORK( J )*SQRT( TEMP ) END IF END IF 30 CONTINUE * 40 CONTINUE END IF RETURN * * End of CGEQPF * END ********************************************************************* SUBROUTINE CGEQRF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) * * -- LAPACK routine (version 1.0b) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. INTEGER INFO, LDA, LWORK, M, N * .. * .. Array Arguments .. COMPLEX A( LDA, * ), TAU( * ), WORK( LWORK ) * .. * * Purpose * ======= * * CGEQRF computes a QR factorization of a complex m by n matrix A: * A = Q * R. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) COMPLEX array, dimension (LDA,N) * On entry, the m by n matrix A. * On exit, the elements on and above the diagonal of the array * contain the min(m,n) by n upper trapezoidal matrix R (R is * upper triangular if m >= n); the elements below the diagonal, * with the array TAU, represent the unitary matrix Q as a * product of elementary reflectors (see Further Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) COMPLEX array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * WORK (workspace) COMPLEX array, dimension (LWORK) * On exit, if INFO = 0, WORK(1) returns the minimum value of * LWORK required to use the optimal blocksize. * * LWORK (input) INTEGER * The dimension of the array WORK. LWORK >= max(1,N). * For optimum performance LWORK should be at least N*NB, * where NB is the optimal blocksize. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * Further Details * =============== * * The matrix Q is represented as a product of elementary reflectors * * Q = H(1) H(2) . . . H(k), where k = min(m,n). * * Each H(i) has the form * * H(i) = I - tau * v * v' * * where tau is a complex scalar, and v is a complex vector with * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), * and tau in TAU(i). * * ===================================================================== * * .. Local Scalars .. INTEGER I, IB, IINFO, IWS, K, LDWORK, NB, NBMIN, NX * .. * .. External Subroutines .. EXTERNAL CGEQR2, CLARFB, CLARFT, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. External Functions .. INTEGER ILAENV EXTERNAL ILAENV * .. * .. Executable Statements .. * * Test the input arguments * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 ELSE IF( LWORK.LT.MAX( 1, N ) ) THEN INFO = -7 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CGEQRF', -INFO ) RETURN END IF * * Quick return if possible * K = MIN( M, N ) IF( K.EQ.0 ) THEN WORK( 1 ) = 1 RETURN END IF * * Determine the block size. * NB = ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) NBMIN = 2 NX = 0 IWS = N IF( NB.GT.1 .AND. NB.LT.K ) THEN * * Determine when to cross over from blocked to unblocked code. * NX = MAX( 0, ILAENV( 3, 'CGEQRF', ' ', M, N, -1, -1 ) ) IF( NX.LT.K ) THEN * * Determine if workspace is large enough for blocked code. * LDWORK = N IWS = LDWORK*NB IF( LWORK.LT.IWS ) THEN * * Not enough workspace to use optimal NB: reduce NB and * determine the minimum value of NB. * NB = LWORK / LDWORK NBMIN = MAX( 2, ILAENV( 2, 'CGEQRF', ' ', M, N, -1, $ -1 ) ) END IF END IF END IF * IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN * * Use blocked code initially * DO 10 I = 1, K - NX, NB IB = MIN( K-I+1, NB ) * * Compute the QR factorization of the current block * A(i:m,i:i+ib-1) * CALL CGEQR2( M-I+1, IB, A( I, I ), LDA, TAU( I ), WORK, $ IINFO ) IF( I+IB.LE.N ) THEN * * Form the triangular factor of the block reflector * H = H(i) H(i+1) . . . H(i+ib-1) * CALL CLARFT( 'Forward', 'Columnwise', M-I+1, IB, $ A( I, I ), LDA, TAU( I ), WORK, LDWORK ) * * Apply H' to A(i:m,i+ib:n) from the left * CALL CLARFB( 'Left', 'Conjugate transpose', 'Forward', $ 'Columnwise', M-I+1, N-I-IB+1, IB, $ A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ), $ LDA, WORK( IB+1 ), LDWORK ) END IF 10 CONTINUE ELSE I = 1 END IF * * Use unblocked code to factor the last or only block. * IF( I.LE.K ) $ CALL CGEQR2( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK, $ IINFO ) * WORK( 1 ) = IWS RETURN * * End of CGEQRF * END ********************************************************************* SUBROUTINE CGEQR2( M, N, A, LDA, TAU, WORK, INFO ) * * -- LAPACK routine (version 1.0b) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. COMPLEX A( LDA, * ), TAU( * ), WORK( * ) * .. * * Purpose * ======= * * CGEQR2 computes a QR factorization of a complex m by n matrix A: * A = Q * R. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) COMPLEX array, dimension (LDA,N) * On entry, the m by n matrix A. * On exit, the elements on and above the diagonal of the array * contain the min(m,n) by n upper trapezoidal matrix R (R is * upper triangular if m >= n); the elements below the diagonal, * with the array TAU, represent the unitary matrix Q as a * product of elementary reflectors (see Further Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) COMPLEX array, dimension (min(M,N)) * The scalar factors of the elementary reflectors (see Further * Details). * * WORK (workspace) COMPLEX array, dimension (N) * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * Further Details * =============== * * The matrix Q is represented as a product of elementary reflectors * * Q = H(1) H(2) . . . H(k), where k = min(m,n). * * Each H(i) has the form * * H(i) = I - tau * v * v' * * where tau is a complex scalar, and v is a complex vector with * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), * and tau in TAU(i). * * ===================================================================== * * .. Parameters .. COMPLEX ONE PARAMETER ( ONE = 1.0E+0 ) * .. * .. Local Scalars .. INTEGER I, K COMPLEX ALPHA * .. * .. External Subroutines .. EXTERNAL CLARF, CLARFG, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC CONJG, MAX, MIN * .. * .. Executable Statements .. * * Test the input arguments * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CGEQR2', -INFO ) RETURN END IF * K = MIN( M, N ) * DO 10 I = 1, K * * Generate elementary reflector H(i) to annihilate A(i+1:m,i) * CALL CLARFG( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1, $ TAU( I ) ) IF( I.LT.N ) THEN * * Apply H(i)' to A(i:m,i+1:n) from the left * ALPHA = A( I, I ) A( I, I ) = ONE CALL CLARF( 'Left', M-I+1, N-I, A( I, I ), 1, $ CONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK ) A( I, I ) = ALPHA END IF 10 CONTINUE RETURN * * End of CGEQR2 * END ********************************************************************* SHAR_EOF fi # end of overwriting check if test -f 'ctrqpx.f' then echo shar: will not over-write existing file "'ctrqpx.f'" else cat << SHAR_EOF > 'ctrqpx.f' SUBROUTINE CTRQPX( JOB, M, N, K, A, LDA, C, LDC, JPVT, IRCOND, $ ORCOND, RANK, SVLUES, WORK, RWORK, INFO ) * * This code is part of a release of the package for computing * rank-revealing QR Factorizations written by: * ================================================================== * Christian H. Bischof and Gregorio Quintana-Orti * Math. and Comp. Sci. Div. Departamento de Informatica * Argonne National Lab. Universidad Jaime I * Argonne, IL 60439 Campus P. Roja, 12071 Castellon * USA Spain * bischof@mcs.anl.gov gquintan@inf.uji.es * ================================================================== * $Revision: 1.42 $ * $Date: 96/12/30 16:59:42 $ * * .. Scalar Arguments .. INTEGER JOB, M, N, K, LDA, LDC, RANK, INFO REAL IRCOND, ORCOND * .. * .. Array Arguments .. COMPLEX A( LDA, * ), C( LDC, * ), WORK( * ) REAL SVLUES( 4 ), RWORK( * ) INTEGER JPVT( * ) * .. * * Purpose * ======= * * CTRQPX detects the right rank for upper triangular matrix A. * The algorithm used here is related to Chandrasekaran&Ipsen * algorithm Hybrid-III. * This algorithm is applied to matrix A until the right rank is * obtained. If the input ordering of matrix A is not accepted, the * matrix will be permuted and retriangularized until the rank is * revealed. * * Arguments * ========= * * JOB (input) INTEGER * The job to do: * = 1: The transformations needed in the * triangularization are only applied to matrix A. * Thus, matrix C is not updated. * = 2: The same transformations needed in the * triangularization of matrix A are applied to * matrix C from the left. * That is, if Q'*A*P=R, then C := Q'*C. * In this case, matrix C is m-by-k. * = 3: The transpose of the transformations needed * in the triangularization of matrix A are applied * to matrix C from the right. * That is, if Q'*A*P=R, then C := C*Q. * In this case, matrix C is k-by-m. * In these three cases, the permutations are always stored * in vector JPVT. * * M (input) INTEGER * The number of rows of matrices A. M >= 0. * If JOB=2, M is the number of rows of matrix C. * If JOB=3, M is the number of columns of matrix C. * * N (input) INTEGER * The number of columns of matrix A. N >= 0. * * K (input) INTEGER * It defines the dimension of matrix C. K >= 0. * If JOB=2, K is the number of columns of matrix C. * If JOB=3, K is the number of rows of matrix C. * * A (input/output) COMPLEX array, dimension (LDA,N) * On entry, the m by n matrix A. * On exit, the upper triangle of the array contains the * min(m,n) by n upper trapezoidal matrix R; the lower triangle * array is filled with zeros. * * LDA (input) INTEGER * The leading dimension of array A. LDA >= max(1,M). * * C (input/output) COMPLEX array, dimension * ( LDC, K ) if JOB=2. * ( LDC, M ) if JOB=3. * If argument JOB asks, all the transformations * applied to matrix A are also applied to matrix C. * * LDC (input) INTEGER * The leading dimension of array C. * If JOB=2, then LDC >= MAX(1,M). * If JOB=3, then LDC >= MAX(1,K). * * JPVT (input/output) INTEGER array, dimension ( N ) * If JPVT(I) = K, then the Ith column of the permuted * A was the Kth column of the original A (just before * the preprocessing). If a permutation occurs, JPVT will * be updated correctly. * JPVT(1:RANK) contains the indices of the columns considered * linearly independent. * JPVT(RANK+1:N) contains the indices of the columns considered * linearly dependent from the previous ones. * * IRCOND (input) FLOATING_DECLARE * 1/IRCOND specifies an upper bound on the condition number * of R11. If IRCOND == 0, IRCOND = machine precision is chosen * as default. IRCOND must be >= 0. * * ORCOND (output) FLOATING_DECLARE * 1/ORCOND is an estimate for the condition number of R11. * * RANK (output) INTEGER * An estimate of the rank offered by this algorithm. * 0 <= RANK <= MIN(M,N). * * SVLUES (output) REAL array, dimension (4) * On exit, SVLUES contains estimates of some of the singular * values of the triangular factor R. * SVLUES(1): largest singular value of R(1:RANK,1:RANK) * SVLUES(2): smallest singular value of R(1:RANK,1:RANK) * SVLUES(3): smallest singular value of R(1:RANK+1,1:RANK+1) * SVLUES(4): smallest singular value of R * If the triangular factorization is a rank-revealing one * (which will be the case if the leading columns were well- * conditioned), then SVLUES(1) will also be an estimate * for the largest singular value of A, SVLUES(2) and SVLUES(3) * will be estimates for the RANK-th and (RANK+1)-st singular * value of A, and SVLUES(4) will be an estimate for the * smallest singular value of A. * By examining these values, one can confirm that the rank is * well defined with respect to the threshold chosen. * * WORK (workspace) COMPLEX array, dimension ( 2*MIN(M,N) ) * * RWORK (workspace) REAL array, dimension ( MIN(M,N)+N ) * * INFO (output) INTEGER * = 0: Successful exit. * < 0: If INFO = -i, the i-th argument had an illegal value * > 0: Problems in the computation of the rank. * 1: Exceeded the allowed maximum number of steps. * 2: Rank not well defined. * In adition, vector SVLUES tell if rank is not well defined. * When INFO.NE.0, the contents of ORCOND may be not the right * one. * * * =================================================================== * * .. Parameters .. INTEGER INB REAL ZERO PARAMETER ( INB = 1, ZERO = 0.0E+0 ) * * Indices into the 'svlues' array. * INTEGER IMAX, IBEFOR, IAFTER, IMIN PARAMETER ( IMAX = 1, IBEFOR = 2, IAFTER = 3, IMIN = 4 ) * .. * .. * .. Common Block .. INTEGER NB COMMON /BSPRQR/ NB * .. * .. Local Scalars .. REAL RCNR, RCNRP1, RCOND LOGICAL GOLEFT, RNKDTD INTEGER MN, OINFO * .. * .. External Functions .. INTEGER ILAENV REAL SLAMCH EXTERNAL ILAENV, SLAMCH * .. * .. External Subroutines .. EXTERNAL XERBLA, CTRQXC * .. * .. Intrinsic Functions .. INTRINSIC MIN, MAX * .. * .. Executable Statements .. * MN = MIN( M, N ) NB = ILAENV( INB, 'CGEQRF', ' ', M, N, 0, 0 ) * * Test input arguments * ==================== * INFO = 0 IF( ( JOB.LT.1 ).OR.( JOB.GT.3 ) ) THEN INFO = -1 ELSE IF( M.LT.0 ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( K.LT.0 ) THEN INFO = -4 ELSE IF( LDA.LT.MAX(1,M) ) THEN INFO = -6 ELSE IF( ( ( JOB.EQ.1 ).AND.( LDC.LT.1 ) ).OR. $ ( ( JOB.EQ.2 ).AND.( LDC.LT.MAX( 1, M ) ) ).OR. $ ( ( JOB.EQ.3 ).AND.( LDC.LT.MAX( 1, K ) ) ) ) THEN INFO = -8 ELSE IF( IRCOND.LT.ZERO ) THEN INFO = -10 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CTRQPX', -INFO ) RETURN END IF * * Quick return if possible. * IF( MN.EQ.0 ) THEN RANK = 0 ORCOND = ZERO SVLUES( IMAX ) = ZERO SVLUES( IBEFOR ) = ZERO SVLUES( IAFTER ) = ZERO SVLUES( IMIN ) = ZERO RETURN END IF * * Check whether Threshold for condition number was supplied. * If not, choose machine precision as default for RCOND. * IF( IRCOND.GT.ZERO ) THEN RCOND = IRCOND ELSE RCOND = SLAMCH( 'Epsilon' ) END IF * * Compute the initial estimate for the rank. * CALL CTRRNK( MN, A, LDA, RCOND, RANK, WORK, INFO ) * * ************************ * * First call to xTRQXC * * ************************ * * Get tighter bounds for the value RANK. * CALL CTRQXC( JOB, M, N, K, A, LDA, C, LDC, JPVT, $ RANK, SVLUES, RCNR, RCNRP1, WORK, RWORK, INFO ) OINFO = 0 IF( INFO.NE.0 ) $ OINFO = INFO * * Check if the numerical rank is larger, equal or smaller than * the contents of RANK. * IF( ( ( RCNR.GE.RCOND ).AND.( RANK.EQ.MN ) ).OR. $ ( ( RCNR.GE.RCOND ).AND.( RCNRP1.LT.RCOND ) ) ) THEN RNKDTD = .TRUE. ELSE IF( ( RCNR.GE.RCOND ).AND.( RCNRP1.GE.RCOND ) ) THEN RNKDTD = .FALSE. GOLEFT = .FALSE. RANK = RANK + 1 ELSE IF( ( RCNR.LT.RCOND ).AND.( RCNRP1.LT.RCOND ) ) THEN IF( RANK.EQ.1 ) THEN RNKDTD = .TRUE. IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN RANK = 0 ORCOND = ZERO SVLUES( IMAX ) = ZERO SVLUES( IBEFOR ) = ZERO SVLUES( IAFTER ) = ZERO SVLUES( IMIN ) = ZERO ELSE RANK = 1 END IF ELSE RNKDTD = .FALSE. GOLEFT = .TRUE. RANK = RANK - 1 END IF ELSE RNKDTD = .TRUE. INFO = 2 END IF * * ***************** * * Start of Loop * * ***************** * * Loop for the detection of the actual rank. The variable RANK is * updated until the rank is found. To avoid infinite loops, the * variable RANK either increases or decreases. * 10 CONTINUE IF( .NOT. RNKDTD ) THEN * * Call to xTRQXC to get tighter bounds for the value RANK. * CALL CTRQXC( JOB, M, N, K, A, LDA, C, LDC, JPVT, $ RANK, SVLUES, RCNR, RCNRP1, WORK, RWORK, INFO ) IF( INFO.NE.0 ) $ OINFO = INFO * * Check if the numerical rank is larger, equal or smaller than * the contents of RANK. * IF( ( ( RCNR.GE.RCOND ).AND.( RANK.EQ.MN ) ).OR. $ ( ( RCNR.GE.RCOND ).AND.( RCNRP1.LT.RCOND ) ) ) THEN RNKDTD = .TRUE. ELSE IF( ( RCNR.GE.RCOND ).AND.( RCNRP1.GE.RCOND ) ) THEN IF( .NOT. GOLEFT ) THEN RANK = RANK + 1 ELSE RNKDTD = .TRUE. INFO = 2 END IF ELSE IF( ( RCNR.LT.RCOND ).AND.( RCNRP1.LT.RCOND ) ) THEN IF( RANK.EQ.1 ) THEN RNKDTD = .TRUE. IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN RANK = 0 ORCOND = ZERO SVLUES( IMAX ) = ZERO SVLUES( IBEFOR ) = ZERO SVLUES( IAFTER ) = ZERO SVLUES( IMIN ) = ZERO ELSE RANK = 1 END IF ELSE GOLEFT = .TRUE. RANK = RANK - 1 END IF ELSE RNKDTD = .TRUE. INFO = 2 END IF * * Jump to the beginning of the loop. * GOTO 10 END IF * * *************** * * end of loop * * *************** * * Give back the obtained value of RCOND and check the value of INFO. * ORCOND = RCNR IF( OINFO.NE.0 ) $ INFO = OINFO * RETURN * * End of CTRQPX * END SHAR_EOF fi # end of overwriting check if test -f 'ctrqpy.f' then echo shar: will not over-write existing file "'ctrqpy.f'" else cat << SHAR_EOF > 'ctrqpy.f' SUBROUTINE CTRQPY( JOB, M, N, K, A, LDA, C, LDC, JPVT, IRCOND, $ ORCOND, RANK, SVLUES, WORK, RWORK, INFO ) * * This code is part of a release of the package for computing * rank-revealing QR Factorizations written by: * ================================================================== * Christian H. Bischof and Gregorio Quintana-Orti * Math. and Comp. Sci. Div. Departamento de Informatica * Argonne National Lab. Universidad Jaime I * Argonne, IL 60439 Campus P. Roja, 12071 Castellon * USA Spain * bischof@mcs.anl.gov gquintan@inf.uji.es * ================================================================== * $Revision: 1.42 $ * $Date: 96/12/30 16:59:43 $ * * .. Scalar Arguments .. INTEGER JOB, M, N, K, LDA, LDC, RANK, INFO REAL IRCOND, ORCOND * .. * .. Array Arguments .. COMPLEX A( LDA, * ), C( LDC, * ), WORK( * ) REAL SVLUES( 4 ), RWORK( * ) INTEGER JPVT( * ) * .. * * Purpose * ======= * * CTRQPY detects the right rank for upper triangular matrix A. * The algorithm used here is an version of Pan and Tang's RRQR * algorithm number 3. * This algorithm is applied to matrix A until the right rank is * obtained. If the input ordering of matrix A is not accepted, the * matrix will be permuted and retriangularized until the rank is * revealed. * * Arguments * ========= * * JOB (input) INTEGER * The job to do: * = 1: The transformations needed in the * triangularization are only applied to matrix A. * Thus, matrix C is not updated. * = 2: The same transformations needed in the * triangularization of matrix A are applied to * matrix C from the left. * That is, if Q'*A*P=R, then C := Q'*C. * In this case, matrix C is m-by-k. * = 3: The transpose of the transformations needed * in the triangularization of matrix A are applied * to matrix C from the right. * That is, if Q'*A*P=R, then C := C*Q. * In this case, matrix C is k-by-m. * In these three cases, the permutations are always stored * in vector JPVT. * * M (input) INTEGER * The number of rows of matrices A. M >= 0. * If JOB=2, M is the number of rows of matrix C. * If JOB=3, M is the number of columns of matrix C. * * N (input) INTEGER * The number of columns of matrix A. N >= 0. * * K (input) INTEGER * It defines the dimension of matrix C. K >= 0. * If JOB=2, K is the number of columns of matrix C. * If JOB=3, K is the number of rows of matrix C. * * A (input/output) COMPLEX array, dimension (LDA,N) * On entry, the m by n matrix A. * On exit, the upper triangle of the array contains the * min(m,n) by n upper trapezoidal matrix R; the lower triangle * array is filled with zeros. * * LDA (input) INTEGER * The leading dimension of array A. LDA >= max(1,M). * * C (input/output) COMPLEX array, dimension * ( LDC, K ) if JOB=2. * ( LDC, M ) if JOB=3. * If argument JOB asks, all the transformations * applied to matrix A are also applied to matrix C. * * LDC (input) INTEGER * The leading dimension of array C. * If JOB=2, then LDC >= MAX(1,M). * If JOB=3, then LDC >= MAX(1,K). * * JPVT (input/output) INTEGER array, dimension ( N ) * If JPVT(I) = K, then the Ith column of the permuted * A was the Kth column of the original A (just before * the preprocessing). If a permutation occurs, JPVT will * be updated correctly. * JPVT(1:RANK) contains the indices of the columns considered * linearly independent. * JPVT(RANK+1:N) contains the indices of the columns considered * linearly dependent from the previous ones. * * IRCOND (input) FLOATING_DECLARE * 1/IRCOND specifies an upper bound on the condition number * of R11. If IRCOND == 0, IRCOND = machine precision is chosen * as default. IRCOND must be >= 0. * * ORCOND (output) FLOATING_DECLARE * 1/ORCOND is an estimate for the condition number of R11. * * RANK (output) INTEGER * An estimate of the rank offered by this algorithm. * 0 <= RANK <= MIN(M,N). * * SVLUES (output) REAL array, dimension (4) * On exit, SVLUES contains estimates of some of the singular * values of the triangular factor R. * SVLUES(1): largest singular value of R(1:RANK,1:RANK) * SVLUES(2): smallest singular value of R(1:RANK,1:RANK) * SVLUES(3): smallest singular value of R(1:RANK+1,1:RANK+1) * SVLUES(4): smallest singular value of R * If the triangular factorization is a rank-revealing one * (which will be the case if the leading columns were well- * conditioned), then SVLUES(1) will also be an estimate * for the largest singular value of A, SVLUES(2) and SVLUES(3) * will be estimates for the RANK-th and (RANK+1)-st singular * value of A, and SVLUES(4) will be an estimate for the * smallest singular value of A. * By examining these values, one can confirm that the rank is * well defined with respect to the threshold chosen. * * WORK (workspace) COMPLEX array, dimension ( MIN(M,N) ) * * RWORK (workspace) REAL array, dimension ( MIN(M,N)+N ) * * INFO (output) INTEGER * = 0: Successful exit. * < 0: If INFO = -i, the i-th argument had an illegal value * > 0: Problems in the computation of the rank. * 1: Exceeded the allowed maximum number of steps. * 2: Rank not well defined. * In adition, vector SVLUES tell if rank is not well defined. * When INFO.NE.0, the contents of ORCOND may be not the right * one. * * * =================================================================== * * .. Parameters .. INTEGER INB REAL ZERO PARAMETER ( INB = 1, ZERO = 0.0E+0 ) * * Indices into the 'svlues' array. * INTEGER IMAX, IBEFOR, IAFTER, IMIN PARAMETER ( IMAX = 1, IBEFOR = 2, IAFTER = 3, IMIN = 4 ) * .. * .. * .. Common Block .. INTEGER NB COMMON /BSPRQR/ NB * .. * .. Local Scalars .. REAL RCNR, RCNRP1, RCOND LOGICAL GOLEFT, RNKDTD INTEGER MN, OINFO * .. * .. External Functions .. INTEGER ILAENV REAL SLAMCH EXTERNAL ILAENV, SLAMCH * .. * .. External Subroutines .. EXTERNAL XERBLA, CTRQYC, CTRRNK * .. * .. Intrinsic Functions .. INTRINSIC MIN, MAX * .. * .. Executable Statements .. * MN = MIN( M, N ) NB = ILAENV( INB, 'CGEQRF', ' ', M, N, 0, 0 ) * * Test input arguments * ==================== * INFO = 0 IF( ( JOB.LT.1 ).OR.( JOB.GT.3 ) ) THEN INFO = -1 ELSE IF( M.LT.0 ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( K.LT.0 ) THEN INFO = -4 ELSE IF( LDA.LT.MAX(1,M) ) THEN INFO = -6 ELSE IF( ( ( JOB.EQ.1 ).AND.( LDC.LT.1 ) ).OR. $ ( ( JOB.EQ.2 ).AND.( LDC.LT.MAX( 1, M ) ) ).OR. $ ( ( JOB.EQ.3 ).AND.( LDC.LT.MAX( 1, K ) ) ) ) THEN INFO = -8 ELSE IF( IRCOND.LT.ZERO ) THEN INFO = -10 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CTRQPY', -INFO ) RETURN END IF * * Quick return if possible. * IF( MN.EQ.0 ) THEN RANK = 0 ORCOND = ZERO SVLUES( IMAX ) = ZERO SVLUES( IBEFOR ) = ZERO SVLUES( IAFTER ) = ZERO SVLUES( IMIN ) = ZERO RETURN END IF * * Check whether Threshold for condition number was supplied. * If not, choose machine precision as default for RCOND. * IF( IRCOND.GT.ZERO ) THEN RCOND = IRCOND ELSE RCOND = SLAMCH( 'Epsilon' ) END IF * * Compute the initial estimate for the rank. * CALL CTRRNK( MN, A, LDA, RCOND, RANK, WORK, INFO ) * * ************************ * * First call to xTRQYC * * ************************ * * Get tighter bounds for the value RANK. * CALL CTRQYC( JOB, M, N, K, A, LDA, C, LDC, JPVT, $ RANK, SVLUES, RCNR, RCNRP1, WORK, RWORK, INFO ) OINFO = 0 IF( INFO.NE.0 ) $ OINFO = INFO * * Check if the numerical rank is larger, equal or smaller than * the contents of RANK. * IF( ( ( RCNR.GE.RCOND ).AND.( RANK.EQ.MN ) ).OR. $ ( ( RCNR.GE.RCOND ).AND.( RCNRP1.LT.RCOND ) ) ) THEN RNKDTD = .TRUE. ELSE IF( ( RCNR.GE.RCOND ).AND.( RCNRP1.GE.RCOND ) ) THEN RNKDTD = .FALSE. GOLEFT = .FALSE. RANK = RANK + 1 ELSE IF( ( RCNR.LT.RCOND ).AND.( RCNRP1.LT.RCOND ) ) THEN IF( RANK.EQ.1 ) THEN RNKDTD = .TRUE. IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN RANK = 0 ORCOND = ZERO SVLUES( IMAX ) = ZERO SVLUES( IBEFOR ) = ZERO SVLUES( IAFTER ) = ZERO SVLUES( IMIN ) = ZERO ELSE RANK = 1 END IF ELSE RNKDTD = .FALSE. GOLEFT = .TRUE. RANK = RANK - 1 END IF ELSE RNKDTD = .TRUE. INFO = 2 END IF * * ***************** * * Start of Loop * * ***************** * * Loop for the detection of the actual rank. The variable RANK is * updated until the rank is found. To avoid infinite loops, the * variable RANK either increases or decreases. * 10 CONTINUE IF( .NOT. RNKDTD ) THEN * * Call to xTRQYC to get tighter bounds for the value RANK. * CALL CTRQYC( JOB, M, N, K, A, LDA, C, LDC, JPVT, $ RANK, SVLUES, RCNR, RCNRP1, WORK, RWORK, INFO ) IF( INFO.NE.0 ) $ OINFO = INFO * * Check if the numerical rank is larger, equal or smaller than * the contents of RANK. * IF( ( ( RCNR.GE.RCOND ).AND.( RANK.EQ.MN ) ).OR. $ ( ( RCNR.GE.RCOND ).AND.( RCNRP1.LT.RCOND ) ) ) THEN RNKDTD = .TRUE. ELSE IF( ( RCNR.GE.RCOND ).AND.( RCNRP1.GE.RCOND ) ) THEN IF( .NOT. GOLEFT ) THEN RANK = RANK + 1 ELSE RNKDTD = .TRUE. INFO = 2 END IF ELSE IF( ( RCNR.LT.RCOND ).AND.( RCNRP1.LT.RCOND ) ) THEN IF( RANK.EQ.1 ) THEN RNKDTD = .TRUE. IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN RANK = 0 ORCOND = ZERO SVLUES( IMAX ) = ZERO SVLUES( IBEFOR ) = ZERO SVLUES( IAFTER ) = ZERO SVLUES( IMIN ) = ZERO ELSE RANK = 1 END IF ELSE GOLEFT = .TRUE. RANK = RANK - 1 END IF ELSE RNKDTD = .TRUE. INFO = 2 END IF * * Jump to the beginning of the loop. * GOTO 10 END IF * * *************** * * end of loop * * *************** * * Give back the obtained value of RCOND and check the value of INFO. * ORCOND = RCNR IF( OINFO.NE.0 ) $ INFO = OINFO * RETURN * * End of CTRQPY * END SHAR_EOF fi # end of overwriting check if test -f 'ctrqxc.f' then echo shar: will not over-write existing file "'ctrqxc.f'" else cat << SHAR_EOF > 'ctrqxc.f' SUBROUTINE CTRQXC( JOB, M, N, K, A, LDA, C, LDC, JPVT, $ RANK, SVLUES, RCNR, RCNRP1, WORK, RWORK, INFO ) * * This code is part of a release of the package for computing * rank-revealing QR Factorizations written by: * ================================================================== * Christian H. Bischof and Gregorio Quintana-Orti * Math. and Comp. Sci. Div. Departamento de Informatica * Argonne National Lab. Universidad Jaime I * Argonne, IL 60439 Campus P. Roja, 12071 Castellon * USA Spain * bischof@mcs.anl.gov gquintan@inf.uji.es * ================================================================== * $Revision: 1.42 $ * $Date: 96/12/30 16:59:43 $ * * .. Scalar Arguments .. INTEGER JOB, M, N, K, LDA, LDC, RANK, INFO REAL RCNR, RCNRP1 * .. * .. Array Arguments .. COMPLEX A( LDA, * ), C( LDC, * ), WORK( * ) REAL SVLUES( 4 ), RWORK( * ) INTEGER