diff --git a/lapack-netlib/SRC/CMakeLists.txt b/lapack-netlib/SRC/CMakeLists.txt index 944401beb..f19bdd302 100644 --- a/lapack-netlib/SRC/CMakeLists.txt +++ b/lapack-netlib/SRC/CMakeLists.txt @@ -106,7 +106,7 @@ set(SLASRC slatbs.f slatdf.f slatps.f slatrd.f slatrs.f slatrz.f slauu2.f slauum.f sopgtr.f sopmtr.f sorg2l.f sorg2r.f sorgbr.f sorghr.f sorgl2.f sorglq.f sorgql.f sorgqr.f sorgr2.f - sorgrq.f sorgtr.f sorm2l.f sorm2r.f sorm22.f + sorgrq.f sorgtr.f sorgtsqr.f sorm2l.f sorm2r.f sorm22.f sormbr.f sormhr.f sorml2.f sormlq.f sormql.f sormqr.f sormr2.f sormr3.f sormrq.f sormrz.f sormtr.f spbcon.f spbequ.f spbrfs.f spbstf.f spbsv.f spbsvx.f @@ -148,9 +148,11 @@ set(SLASRC sgetsls.f sgeqr.f slatsqr.f slamtsqr.f sgemqr.f sgelq.f slaswlq.f slamswlq.f sgemlq.f stplqt.f stplqt2.f stpmlqt.f + sorhr_col.f slaorhr_col_getrfnp.f slaorhr_col_getrfnp2.f ssytrd_2stage.f ssytrd_sy2sb.f ssytrd_sb2st.F ssb2st_kernels.f ssyevd_2stage.f ssyev_2stage.f ssyevx_2stage.f ssyevr_2stage.f - ssbev_2stage.f ssbevx_2stage.f ssbevd_2stage.f ssygv_2stage.f) + ssbev_2stage.f ssbevx_2stage.f ssbevd_2stage.f ssygv_2stage.f + sgesvdq.f scombssq.f) set(DSLASRC spotrs.f sgetrs.f spotrf.f sgetrf.f) @@ -233,7 +235,7 @@ set(CLASRC ctptrs.f ctrcon.f ctrevc.f ctrevc3.f ctrexc.f ctrrfs.f ctrsen.f ctrsna.f ctrsyl.f ctrti2.f ctrtri.f ctrtrs.f ctzrzf.f cung2l.f cung2r.f cungbr.f cunghr.f cungl2.f cunglq.f cungql.f cungqr.f cungr2.f - cungrq.f cungtr.f cunm2l.f cunm2r.f cunmbr.f cunmhr.f cunml2.f cunm22.f + cungrq.f cungtr.f cungtsqr.f cunm2l.f cunm2r.f cunmbr.f cunmhr.f cunml2.f cunm22.f cunmlq.f cunmql.f cunmqr.f cunmr2.f cunmr3.f cunmrq.f cunmrz.f cunmtr.f cupgtr.f cupmtr.f icmax1.f scsum1.f cstemr.f chfrk.f ctfttp.f clanhf.f cpftrf.f cpftri.f cpftrs.f ctfsm.f ctftri.f @@ -247,9 +249,11 @@ set(CLASRC cgetsls.f cgeqr.f clatsqr.f clamtsqr.f cgemqr.f cgelq.f claswlq.f clamswlq.f cgemlq.f ctplqt.f ctplqt2.f ctpmlqt.f + cunhr_col.f claunhr_col_getrfnp.f claunhr_col_getrfnp2.f chetrd_2stage.f chetrd_he2hb.f chetrd_hb2st.F chb2st_kernels.f cheevd_2stage.f cheev_2stage.f cheevx_2stage.f cheevr_2stage.f - chbev_2stage.f chbevx_2stage.f chbevd_2stage.f chegv_2stage.f) + chbev_2stage.f chbevx_2stage.f chbevd_2stage.f chegv_2stage.f + cgesvdq.f) set(CXLASRC cgesvxx.f cgerfsx.f cla_gerfsx_extended.f cla_geamv.f cla_gercond_c.f cla_gercond_x.f cla_gerpvgrw.f @@ -295,7 +299,7 @@ set(DLASRC dlatbs.f dlatdf.f dlatps.f dlatrd.f dlatrs.f dlatrz.f dlauu2.f dlauum.f dopgtr.f dopmtr.f dorg2l.f dorg2r.f dorgbr.f dorghr.f dorgl2.f dorglq.f dorgql.f dorgqr.f dorgr2.f - dorgrq.f dorgtr.f dorm2l.f dorm2r.f dorm22.f + dorgrq.f dorgtr.f dorgtsqr.f dorm2l.f dorm2r.f dorm22.f dormbr.f dormhr.f dorml2.f dormlq.f dormql.f dormqr.f dormr2.f dormr3.f dormrq.f dormrz.f dormtr.f dpbcon.f dpbequ.f dpbrfs.f dpbstf.f dpbsv.f dpbsvx.f @@ -339,9 +343,11 @@ set(DLASRC dgetsls.f dgeqr.f dlatsqr.f dlamtsqr.f dgemqr.f dgelq.f dlaswlq.f dlamswlq.f dgemlq.f dtplqt.f dtplqt2.f dtpmlqt.f + dorhr_col.f dlaorhr_col_getrfnp.f dlaorhr_col_getrfnp2.f dsytrd_2stage.f dsytrd_sy2sb.f dsytrd_sb2st.F dsb2st_kernels.f dsyevd_2stage.f dsyev_2stage.f dsyevx_2stage.f dsyevr_2stage.f - dsbev_2stage.f dsbevx_2stage.f dsbevd_2stage.f dsygv_2stage.f) + dsbev_2stage.f dsbevx_2stage.f dsbevd_2stage.f dsygv_2stage.f + dgesvdq.f dcombssq.f) set(DXLASRC dgesvxx.f dgerfsx.f dla_gerfsx_extended.f dla_geamv.f dla_gercond.f dla_gerpvgrw.f dsysvxx.f dsyrfsx.f @@ -424,7 +430,7 @@ set(ZLASRC ztptrs.f ztrcon.f ztrevc.f ztrevc3.f ztrexc.f ztrrfs.f ztrsen.f ztrsna.f ztrsyl.f ztrti2.f ztrtri.f ztrtrs.f ztzrzf.f zung2l.f zung2r.f zungbr.f zunghr.f zungl2.f zunglq.f zungql.f zungqr.f zungr2.f - zungrq.f zungtr.f zunm2l.f zunm2r.f zunmbr.f zunmhr.f zunml2.f zunm22.f + zungrq.f zungtr.f zungtsqr.f zunm2l.f zunm2r.f zunmbr.f zunmhr.f zunml2.f zunm22.f zunmlq.f zunmql.f zunmqr.f zunmr2.f zunmr3.f zunmrq.f zunmrz.f zunmtr.f zupgtr.f zupmtr.f izmax1.f dzsum1.f zstemr.f @@ -440,9 +446,11 @@ set(ZLASRC zgelqt.f zgelqt3.f zgemlqt.f zgetsls.f zgeqr.f zlatsqr.f zlamtsqr.f zgemqr.f zgelq.f zlaswlq.f zlamswlq.f zgemlq.f + zunhr_col.f zlaunhr_col_getrfnp.f zlaunhr_col_getrfnp2.f zhetrd_2stage.f zhetrd_he2hb.f zhetrd_hb2st.F zhb2st_kernels.f zheevd_2stage.f zheev_2stage.f zheevx_2stage.f zheevr_2stage.f - zhbev_2stage.f zhbevx_2stage.f zhbevd_2stage.f zhegv_2stage.f) + zhbev_2stage.f zhbevx_2stage.f zhbevd_2stage.f zhegv_2stage.f + zgesvdq.f) set(ZXLASRC zgesvxx.f zgerfsx.f zla_gerfsx_extended.f zla_geamv.f zla_gercond_c.f zla_gercond_x.f zla_gerpvgrw.f zsysvxx.f zsyrfsx.f @@ -504,7 +512,7 @@ if(USE_XBLAS) endif() target_link_libraries(lapack PRIVATE ${BLAS_LIBRARIES}) -if (${CMAKE_BUILD_TYPE_UPPER} STREQUAL "COVERAGE") +if(_is_coverage_build) target_link_libraries(lapack PRIVATE gcov) add_coverage(lapack) endif() diff --git a/lapack-netlib/SRC/Makefile b/lapack-netlib/SRC/Makefile index 1c276aff6..9f79e20e9 100644 --- a/lapack-netlib/SRC/Makefile +++ b/lapack-netlib/SRC/Makefile @@ -1,5 +1,3 @@ -include ../make.inc - ####################################################################### # This is the makefile to create a library for LAPACK. # The files are organized as follows: @@ -44,7 +42,7 @@ include ../make.inc # and is created at the next higher directory level. # # To remove the object files after the library is created, enter -# make clean +# make cleanobj # On some systems, you can force the source files to be recompiled by # entering (for example) # make single FRC=FRC @@ -56,6 +54,13 @@ include ../make.inc # ####################################################################### +TOPSRCDIR = .. +include $(TOPSRCDIR)/make.inc + +.SUFFIXES: .F .o +.F.o: + $(FC) $(FFLAGS) -c -o $@ $< + ALLAUX_O = ilaenv.o ilaenv2stage.o ieeeck.o lsamen.o xerbla.o xerbla_array.o \ iparmq.o iparam2stage.o \ ilaprec.o ilatrans.o ilauplo.o iladiag.o chla_transtype.o \ @@ -128,7 +133,7 @@ SLASRC_O = \ slatbs.o slatdf.o slatps.o slatrd.o slatrs.o slatrz.o \ slauu2.o slauum.o sopgtr.o sopmtr.o sorg2l.o sorg2r.o \ sorgbr.o sorghr.o sorgl2.o sorglq.o sorgql.o sorgqr.o sorgr2.o \ - sorgrq.o sorgtr.o sorm2l.o sorm2r.o sorm22.o \ + sorgrq.o sorgtr.o sorgtsqr.o sorm2l.o sorm2r.o sorm22.o \ sormbr.o sormhr.o sorml2.o sormlq.o sormql.o sormqr.o sormr2.o \ sormr3.o sormrq.o sormrz.o sormtr.o spbcon.o spbequ.o spbrfs.o \ spbstf.o spbsv.o spbsvx.o \ @@ -171,9 +176,11 @@ SLASRC_O = \ sgetsls.o sgeqr.o slatsqr.o slamtsqr.o sgemqr.o \ sgelq.o slaswlq.o slamswlq.o sgemlq.o \ stplqt.o stplqt2.o stpmlqt.o \ + sorhr_col.o slaorhr_col_getrfnp.o slaorhr_col_getrfnp2.o \ ssytrd_2stage.o ssytrd_sy2sb.o ssytrd_sb2st.o ssb2st_kernels.o \ ssyevd_2stage.o ssyev_2stage.o ssyevx_2stage.o ssyevr_2stage.o \ - ssbev_2stage.o ssbevx_2stage.o ssbevd_2stage.o ssygv_2stage.o + ssbev_2stage.o ssbevx_2stage.o ssbevd_2stage.o ssygv_2stage.o \ + sgesvdq.o scombssq.o DSLASRC_O = spotrs.o sgetrs.o spotrf.o sgetrf.o @@ -258,7 +265,7 @@ CLASRC_O = \ ctptrs.o ctrcon.o ctrevc.o ctrevc3.o ctrexc.o ctrrfs.o ctrsen.o ctrsna.o \ ctrsyl.o ctrti2.o ctrtri.o ctrtrs.o ctzrzf.o cung2l.o cung2r.o \ cungbr.o cunghr.o cungl2.o cunglq.o cungql.o cungqr.o cungr2.o \ - cungrq.o cungtr.o cunm2l.o cunm2r.o cunmbr.o cunmhr.o cunml2.o cunm22.o \ + cungrq.o cungtr.o cungtsqr.o cunm2l.o cunm2r.o cunmbr.o cunmhr.o cunml2.o cunm22.o \ cunmlq.o cunmql.o cunmqr.o cunmr2.o cunmr3.o cunmrq.o cunmrz.o \ cunmtr.o cupgtr.o cupmtr.o icmax1.o scsum1.o cstemr.o \ chfrk.o ctfttp.o clanhf.o cpftrf.o cpftri.o cpftrs.o ctfsm.o ctftri.o \ @@ -272,9 +279,11 @@ CLASRC_O = \ cgetsls.o cgeqr.o clatsqr.o clamtsqr.o cgemqr.o \ cgelq.o claswlq.o clamswlq.o cgemlq.o \ ctplqt.o ctplqt2.o ctpmlqt.o \ + cunhr_col.o claunhr_col_getrfnp.o claunhr_col_getrfnp2.o \ chetrd_2stage.o chetrd_he2hb.o chetrd_hb2st.o chb2st_kernels.o \ cheevd_2stage.o cheev_2stage.o cheevx_2stage.o cheevr_2stage.o \ - chbev_2stage.o chbevx_2stage.o chbevd_2stage.o chegv_2stage.o + chbev_2stage.o chbevx_2stage.o chbevd_2stage.o chegv_2stage.o \ + cgesvdq.o ifdef USEXBLAS CXLASRC = cgesvxx.o cgerfsx.o cla_gerfsx_extended.o cla_geamv.o \ @@ -324,7 +333,7 @@ DLASRC_O = \ dlatbs.o dlatdf.o dlatps.o dlatrd.o dlatrs.o dlatrz.o dlauu2.o \ dlauum.o dopgtr.o dopmtr.o dorg2l.o dorg2r.o \ dorgbr.o dorghr.o dorgl2.o dorglq.o dorgql.o dorgqr.o dorgr2.o \ - dorgrq.o dorgtr.o dorm2l.o dorm2r.o dorm22.o \ + dorgrq.o dorgtr.o dorgtsqr.o dorm2l.o dorm2r.o dorm22.o \ dormbr.o dormhr.o dorml2.o dormlq.o dormql.o dormqr.o dormr2.o \ dormr3.o dormrq.o dormrz.o dormtr.o dpbcon.o dpbequ.o dpbrfs.o \ dpbstf.o dpbsv.o dpbsvx.o \ @@ -368,9 +377,11 @@ DLASRC_O = \ dgetsls.o dgeqr.o dlatsqr.o dlamtsqr.o dgemqr.o \ dgelq.o dlaswlq.o dlamswlq.o dgemlq.o \ dtplqt.o dtplqt2.o dtpmlqt.o \ + dorhr_col.o dlaorhr_col_getrfnp.o dlaorhr_col_getrfnp2.o \ dsytrd_2stage.o dsytrd_sy2sb.o dsytrd_sb2st.o dsb2st_kernels.o \ dsyevd_2stage.o dsyev_2stage.o dsyevx_2stage.o dsyevr_2stage.o \ - dsbev_2stage.o dsbevx_2stage.o dsbevd_2stage.o dsygv_2stage.o + dsbev_2stage.o dsbevx_2stage.o dsbevd_2stage.o dsygv_2stage.o \ + dgesvdq.o dcombssq.o ifdef USEXBLAS DXLASRC = dgesvxx.o dgerfsx.o dla_gerfsx_extended.o dla_geamv.o \ @@ -456,7 +467,7 @@ ZLASRC_O = \ ztptrs.o ztrcon.o ztrevc.o ztrevc3.o ztrexc.o ztrrfs.o ztrsen.o ztrsna.o \ ztrsyl.o ztrti2.o ztrtri.o ztrtrs.o ztzrzf.o zung2l.o \ zung2r.o zungbr.o zunghr.o zungl2.o zunglq.o zungql.o zungqr.o zungr2.o \ - zungrq.o zungtr.o zunm2l.o zunm2r.o zunmbr.o zunmhr.o zunml2.o zunm22.o \ + zungrq.o zungtr.o zungtsqr.o zunm2l.o zunm2r.o zunmbr.o zunmhr.o zunml2.o zunm22.o \ zunmlq.o zunmql.o zunmqr.o zunmr2.o zunmr3.o zunmrq.o zunmrz.o \ zunmtr.o zupgtr.o \ zupmtr.o izmax1.o dzsum1.o zstemr.o \ @@ -472,9 +483,11 @@ ZLASRC_O = \ zgelqt.o zgelqt3.o zgemlqt.o \ zgetsls.o zgeqr.o zlatsqr.o zlamtsqr.o zgemqr.o \ zgelq.o zlaswlq.o zlamswlq.o zgemlq.o \ + zunhr_col.o zlaunhr_col_getrfnp.o zlaunhr_col_getrfnp2.o \ zhetrd_2stage.o zhetrd_he2hb.o zhetrd_hb2st.o zhb2st_kernels.o \ zheevd_2stage.o zheev_2stage.o zheevx_2stage.o zheevr_2stage.o \ - zhbev_2stage.o zhbevx_2stage.o zhbevd_2stage.o zhegv_2stage.o + zhbev_2stage.o zhbevx_2stage.o zhbevd_2stage.o zhegv_2stage.o \ + zgesvdq.o ifdef USEXBLAS ZXLASRC = zgesvxx.o zgerfsx.o zla_gerfsx_extended.o zla_geamv.o \ @@ -550,33 +563,29 @@ ifdef BUILD_DEPRECATED DEPRECATED = $(DEPRECSRC) endif -all: ../$(LAPACKLIB) +.PHONY: all +all: $(LAPACKLIB) -.PHONY: ../$(LAPACKLIB) - -../$(LAPACKLIB): $(ALLOBJ) $(ALLXOBJ) $(DEPRECATED) - $(ARCH) $(ARCHFLAGS) $@ $(ALLOBJ) $(ALLXOBJ) $(DEPRECATED) +$(LAPACKLIB): $(ALLOBJ) $(ALLXOBJ) $(DEPRECATED) + $(AR) $(ARFLAGS) $@ $^ $(RANLIB) $@ +.PHONY: single complex double complex16 single: $(SLASRC) $(DSLASRC) $(SXLASRC) $(SCLAUX) $(ALLAUX) - $(ARCH) $(ARCHFLAGS) ../$(LAPACKLIB) $(SLASRC) $(DSLASRC) \ - $(SXLASRC) $(SCLAUX) $(ALLAUX) - $(RANLIB) ../$(LAPACKLIB) + $(AR) $(ARFLAGS) $(LAPACKLIB) $^ + $(RANLIB) $(LAPACKLIB) complex: $(CLASRC) $(ZCLASRC) $(CXLASRC) $(SCLAUX) $(ALLAUX) - $(ARCH) $(ARCHFLAGS) ../$(LAPACKLIB) $(CLASRC) $(ZCLASRC) \ - $(CXLASRC) $(SCLAUX) $(ALLAUX) - $(RANLIB) ../$(LAPACKLIB) + $(AR) $(ARFLAGS) $(LAPACKLIB) $^ + $(RANLIB) $(LAPACKLIB) double: $(DLASRC) $(DSLASRC) $(DXLASRC) $(DZLAUX) $(ALLAUX) - $(ARCH) $(ARCHFLAGS) ../$(LAPACKLIB) $(DLASRC) $(DSLASRC) \ - $(DXLASRC) $(DZLAUX) $(ALLAUX) - $(RANLIB) ../$(LAPACKLIB) + $(AR) $(ARFLAGS) $(LAPACKLIB) $^ + $(RANLIB) $(LAPACKLIB) complex16: $(ZLASRC) $(ZCLASRC) $(ZXLASRC) $(DZLAUX) $(ALLAUX) - $(ARCH) $(ARCHFLAGS) ../$(LAPACKLIB) $(ZLASRC) $(ZCLASRC) \ - $(ZXLASRC) $(DZLAUX) $(ALLAUX) - $(RANLIB) ../$(LAPACKLIB) + $(AR) $(ARFLAGS) $(LAPACKLIB) $^ + $(RANLIB) $(LAPACKLIB) $(ALLAUX): $(FRC) $(SCLAUX): $(FRC) @@ -597,18 +606,16 @@ endif FRC: @FRC=$(FRC) -clean: +.PHONY: clean cleanobj cleanlib +clean: cleanobj cleanlib +cleanobj: rm -f *.o DEPRECATED/*.o +cleanlib: + rm -f $(LAPACKLIB) -.f.o: - $(FORTRAN) $(OPTS) -c -o $@ $< - -.F.o: - $(FORTRAN) $(OPTS) -c $< -o $@ - -slaruv.o: slaruv.f ; $(FORTRAN) $(NOOPT) -c -o $@ $< -dlaruv.o: dlaruv.f ; $(FORTRAN) $(NOOPT) -c -o $@ $< -sla_wwaddw.o: sla_wwaddw.f ; $(FORTRAN) $(NOOPT) -c -o $@ $< -dla_wwaddw.o: dla_wwaddw.f ; $(FORTRAN) $(NOOPT) -c -o $@ $< -cla_wwaddw.o: cla_wwaddw.f ; $(FORTRAN) $(NOOPT) -c -o $@ $< -zla_wwaddw.o: zla_wwaddw.f ; $(FORTRAN) $(NOOPT) -c -o $@ $< +slaruv.o: slaruv.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< +dlaruv.o: dlaruv.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< +sla_wwaddw.o: sla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< +dla_wwaddw.o: dla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< +cla_wwaddw.o: cla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< +zla_wwaddw.o: zla_wwaddw.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< diff --git a/lapack-netlib/SRC/cgbrfsx.f b/lapack-netlib/SRC/cgbrfsx.f index 041b6a1b6..c23608afb 100644 --- a/lapack-netlib/SRC/cgbrfsx.f +++ b/lapack-netlib/SRC/cgbrfsx.f @@ -75,7 +75,7 @@ *> Specifies the form of the system of equations: *> = 'N': A * X = B (No transpose) *> = 'T': A**T * X = B (Transpose) -*> = 'C': A**H * X = B (Conjugate transpose = Transpose) +*> = 'C': A**H * X = B (Conjugate transpose) *> \endverbatim *> *> \param[in] EQUED @@ -308,7 +308,7 @@ *> information as described below. There currently are up to three *> pieces of information returned for each right-hand side. If *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then -*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most +*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most *> the first (:,N_ERR_BNDS) entries are returned. *> *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith @@ -344,14 +344,14 @@ *> \param[in] NPARAMS *> \verbatim *> NPARAMS is INTEGER -*> Specifies the number of parameters set in PARAMS. If .LE. 0, the +*> Specifies the number of parameters set in PARAMS. If <= 0, the *> PARAMS array is never referenced and default values are used. *> \endverbatim *> *> \param[in,out] PARAMS *> \verbatim *> PARAMS is REAL array, dimension NPARAMS -*> Specifies algorithm parameters. If an entry is .LT. 0.0, then +*> Specifies algorithm parameters. If an entry is < 0.0, then *> that entry will be filled with default value used for that *> parameter. Only positions up to NPARAMS are accessed; defaults *> are used for higher-numbered parameters. @@ -359,9 +359,9 @@ *> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative *> refinement or not. *> Default: 1.0 -*> = 0.0 : No refinement is performed, and no error bounds are +*> = 0.0: No refinement is performed, and no error bounds are *> computed. -*> = 1.0 : Use the double-precision refinement algorithm, +*> = 1.0: Use the double-precision refinement algorithm, *> possibly with doubled-single computations if the *> compilation environment does not support DOUBLE *> PRECISION. diff --git a/lapack-netlib/SRC/cgbsvxx.f b/lapack-netlib/SRC/cgbsvxx.f index 2e113f99c..9f2bbbc1c 100644 --- a/lapack-netlib/SRC/cgbsvxx.f +++ b/lapack-netlib/SRC/cgbsvxx.f @@ -431,7 +431,7 @@ *> information as described below. There currently are up to three *> pieces of information returned for each right-hand side. If *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then -*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most +*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most *> the first (:,N_ERR_BNDS) entries are returned. *> *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith @@ -467,14 +467,14 @@ *> \param[in] NPARAMS *> \verbatim *> NPARAMS is INTEGER -*> Specifies the number of parameters set in PARAMS. If .LE. 0, the +*> Specifies the number of parameters set in PARAMS. If <= 0, the *> PARAMS array is never referenced and default values are used. *> \endverbatim *> *> \param[in,out] PARAMS *> \verbatim *> PARAMS is REAL array, dimension NPARAMS -*> Specifies algorithm parameters. If an entry is .LT. 0.0, then +*> Specifies algorithm parameters. If an entry is < 0.0, then *> that entry will be filled with default value used for that *> parameter. Only positions up to NPARAMS are accessed; defaults *> are used for higher-numbered parameters. @@ -482,9 +482,9 @@ *> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative *> refinement or not. *> Default: 1.0 -*> = 0.0 : No refinement is performed, and no error bounds are +*> = 0.0: No refinement is performed, and no error bounds are *> computed. -*> = 1.0 : Use the double-precision refinement algorithm, +*> = 1.0: Use the double-precision refinement algorithm, *> possibly with doubled-single computations if the *> compilation environment does not support DOUBLE *> PRECISION. diff --git a/lapack-netlib/SRC/cgebak.f b/lapack-netlib/SRC/cgebak.f index 63c73bfa7..9b6402622 100644 --- a/lapack-netlib/SRC/cgebak.f +++ b/lapack-netlib/SRC/cgebak.f @@ -48,10 +48,10 @@ *> \verbatim *> JOB is CHARACTER*1 *> Specifies the type of backward transformation required: -*> = 'N', do nothing, return immediately; -*> = 'P', do backward transformation for permutation only; -*> = 'S', do backward transformation for scaling only; -*> = 'B', do backward transformations for both permutation and +*> = 'N': do nothing, return immediately; +*> = 'P': do backward transformation for permutation only; +*> = 'S': do backward transformation for scaling only; +*> = 'B': do backward transformations for both permutation and *> scaling. *> JOB must be the same as the argument JOB supplied to CGEBAL. *> \endverbatim diff --git a/lapack-netlib/SRC/cgeev.f b/lapack-netlib/SRC/cgeev.f index bdd75e4f1..f07d9b755 100644 --- a/lapack-netlib/SRC/cgeev.f +++ b/lapack-netlib/SRC/cgeev.f @@ -157,7 +157,7 @@ *> < 0: if INFO = -i, the i-th argument had an illegal value. *> > 0: if INFO = i, the QR algorithm failed to compute all the *> eigenvalues, and no eigenvectors have been computed; -*> elements and i+1:N of W contain eigenvalues which have +*> elements i+1:N of W contain eigenvalues which have *> converged. *> \endverbatim * diff --git a/lapack-netlib/SRC/cgejsv.f b/lapack-netlib/SRC/cgejsv.f index a7b1c451c..350da4c40 100644 --- a/lapack-netlib/SRC/cgejsv.f +++ b/lapack-netlib/SRC/cgejsv.f @@ -80,13 +80,13 @@ *> desirable, then this option is advisable. The input matrix A *> is preprocessed with QR factorization with FULL (row and *> column) pivoting. -*> = 'G' Computation as with 'F' with an additional estimate of the +*> = 'G': Computation as with 'F' with an additional estimate of the *> condition number of B, where A=B*D. If A has heavily weighted *> rows, then using this condition number gives too pessimistic *> error bound. *> = 'A': Small singular values are not well determined by the data *> and are considered as noisy; the matrix is treated as -*> numerically rank defficient. The error in the computed +*> numerically rank deficient. The error in the computed *> singular values is bounded by f(m,n)*epsilon*||A||. *> The computed SVD A = U * S * V^* restores A up to *> f(m,n)*epsilon*||A||. @@ -117,7 +117,7 @@ *> = 'V': N columns of V are returned in the array V; Jacobi rotations *> are not explicitly accumulated. *> = 'J': N columns of V are returned in the array V, but they are -*> computed as the product of Jacobi rotations, if JOBT .EQ. 'N'. +*> computed as the product of Jacobi rotations, if JOBT = 'N'. *> = 'W': V may be used as workspace of length N*N. See the description *> of V. *> = 'N': V is not computed. @@ -131,7 +131,7 @@ *> specified range. If A .NE. 0 is scaled so that the largest singular *> value of c*A is around SQRT(BIG), BIG=SLAMCH('O'), then JOBR issues *> the licence to kill columns of A whose norm in c*A is less than -*> SQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN, +*> SQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN, *> where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E'). *> = 'N': Do not kill small columns of c*A. This option assumes that *> BLAS and QR factorizations and triangular solvers are @@ -229,7 +229,7 @@ *> If JOBU = 'F', then U contains on exit the M-by-M matrix of *> the left singular vectors, including an ONB *> of the orthogonal complement of the Range(A). -*> If JOBU = 'W' .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N), +*> If JOBU = 'W' .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N), *> then U is used as workspace if the procedure *> replaces A with A^*. In that case, [V] is computed *> in U as left singular vectors of A^* and then @@ -251,7 +251,7 @@ *> V is COMPLEX array, dimension ( LDV, N ) *> If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of *> the right singular vectors; -*> If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N), +*> If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N), *> then V is used as workspace if the pprocedure *> replaces A with A^*. In that case, [U] is computed *> in V as right singular vectors of A^* and then @@ -282,7 +282,7 @@ *> Length of CWORK to confirm proper allocation of workspace. *> LWORK depends on the job: *> -*> 1. If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and +*> 1. If only SIGMA is needed ( JOBU = 'N', JOBV = 'N' ) and *> 1.1 .. no scaled condition estimate required (JOBA.NE.'E'.AND.JOBA.NE.'G'): *> LWORK >= 2*N+1. This is the minimal requirement. *> ->> For optimal performance (blocked code) the optimal value @@ -298,9 +298,9 @@ *> In general, the optimal length LWORK is computed as *> LWORK >= max(N+LWORK(CGEQP3),N+LWORK(CGEQRF), LWORK(CGESVJ), *> N*N+LWORK(CPOCON)). -*> 2. If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'), -*> (JOBU.EQ.'N') -*> 2.1 .. no scaled condition estimate requested (JOBE.EQ.'N'): +*> 2. If SIGMA and the right singular vectors are needed (JOBV = 'V'), +*> (JOBU = 'N') +*> 2.1 .. no scaled condition estimate requested (JOBE = 'N'): *> -> the minimal requirement is LWORK >= 3*N. *> -> For optimal performance, *> LWORK >= max(N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, @@ -318,10 +318,10 @@ *> LWORK >= max(N+LWORK(CGEQP3), LWORK(CPOCON), N+LWORK(CGESVJ), *> N+LWORK(CGELQF), 2*N+LWORK(CGEQRF), N+LWORK(CUNMLQ)). *> 3. If SIGMA and the left singular vectors are needed -*> 3.1 .. no scaled condition estimate requested (JOBE.EQ.'N'): +*> 3.1 .. no scaled condition estimate requested (JOBE = 'N'): *> -> the minimal requirement is LWORK >= 3*N. *> -> For optimal performance: -*> if JOBU.EQ.'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, +*> if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, *> where NB is the optimal block size for CGEQP3, CGEQRF, CUNMQR. *> In general, the optimal length LWORK is computed as *> LWORK >= max(N+LWORK(CGEQP3), 2*N+LWORK(CGEQRF), N+LWORK(CUNMQR)). @@ -329,16 +329,16 @@ *> required (JOBA='E', or 'G'). *> -> the minimal requirement is LWORK >= 3*N. *> -> For optimal performance: -*> if JOBU.EQ.'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, +*> if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, *> where NB is the optimal block size for CGEQP3, CGEQRF, CUNMQR. *> In general, the optimal length LWORK is computed as *> LWORK >= max(N+LWORK(CGEQP3),N+LWORK(CPOCON), *> 2*N+LWORK(CGEQRF), N+LWORK(CUNMQR)). *> -*> 4. If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and -*> 4.1. if JOBV.EQ.'V' +*> 4. If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and +*> 4.1. if JOBV = 'V' *> the minimal requirement is LWORK >= 5*N+2*N*N. -*> 4.2. if JOBV.EQ.'J' the minimal requirement is +*> 4.2. if JOBV = 'J' the minimal requirement is *> LWORK >= 4*N+N*N. *> In both cases, the allocated CWORK can accommodate blocked runs *> of CGEQP3, CGEQRF, CGELQF, CUNMQR, CUNMLQ. @@ -357,7 +357,7 @@ *> of A. (See the description of SVA().) *> RWORK(2) = See the description of RWORK(1). *> RWORK(3) = SCONDA is an estimate for the condition number of -*> column equilibrated A. (If JOBA .EQ. 'E' or 'G') +*> column equilibrated A. (If JOBA = 'E' or 'G') *> SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1). *> It is computed using SPOCON. It holds *> N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA @@ -376,7 +376,7 @@ *> triangular factor in the first QR factorization. *> RWORK(5) = an estimate of the scaled condition number of the *> triangular factor in the second QR factorization. -*> The following two parameters are computed if JOBT .EQ. 'T'. +*> The following two parameters are computed if JOBT = 'T'. *> They are provided for a developer/implementer who is familiar *> with the details of the method. *> RWORK(6) = the entropy of A^* * A :: this is the Shannon entropy @@ -457,23 +457,23 @@ *> of JOBA and JOBR. *> IWORK(2) = the number of the computed nonzero singular values *> IWORK(3) = if nonzero, a warning message: -*> If IWORK(3).EQ.1 then some of the column norms of A +*> If IWORK(3) = 1 then some of the column norms of A *> were denormalized floats. The requested high accuracy *> is not warranted by the data. -*> IWORK(4) = 1 or -1. If IWORK(4) .EQ. 1, then the procedure used A^* to +*> IWORK(4) = 1 or -1. If IWORK(4) = 1, then the procedure used A^* to *> do the job as specified by the JOB parameters. -*> If the call to CGEJSV is a workspace query (indicated by LWORK .EQ. -1 and -*> LRWORK .EQ. -1), then on exit IWORK(1) contains the required length of +*> If the call to CGEJSV is a workspace query (indicated by LWORK = -1 and +*> LRWORK = -1), then on exit IWORK(1) contains the required length of *> IWORK for the job parameters used in the call. *> \endverbatim *> *> \param[out] INFO *> \verbatim *> INFO is INTEGER -*> < 0 : if INFO = -i, then the i-th argument had an illegal value. -*> = 0 : successful exit; -*> > 0 : CGEJSV did not converge in the maximal allowed number -*> of sweeps. The computed values may be inaccurate. +*> < 0: if INFO = -i, then the i-th argument had an illegal value. +*> = 0: successful exit; +*> > 0: CGEJSV did not converge in the maximal allowed number +*> of sweeps. The computed values may be inaccurate. *> \endverbatim * * Authors: @@ -1336,7 +1336,7 @@ IF ( L2ABER ) THEN * Standard absolute error bound suffices. All sigma_i with * sigma_i < N*EPSLN*||A|| are flushed to zero. This is an -* agressive enforcement of lower numerical rank by introducing a +* aggressive enforcement of lower numerical rank by introducing a * backward error of the order of N*EPSLN*||A||. TEMP1 = SQRT(REAL(N))*EPSLN DO 3001 p = 2, N @@ -1348,9 +1348,9 @@ 3001 CONTINUE 3002 CONTINUE ELSE IF ( L2RANK ) THEN -* .. similarly as above, only slightly more gentle (less agressive). +* .. similarly as above, only slightly more gentle (less aggressive). * Sudden drop on the diagonal of R1 is used as the criterion for -* close-to-rank-defficient. +* close-to-rank-deficient. TEMP1 = SQRT(SFMIN) DO 3401 p = 2, N IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR. @@ -1718,7 +1718,7 @@ CALL CPOCON('L',NR,CWORK(2*N+1),NR,ONE,TEMP1, $ CWORK(2*N+NR*NR+1),RWORK,IERR) CONDR1 = ONE / SQRT(TEMP1) -* .. here need a second oppinion on the condition number +* .. here need a second opinion on the condition number * .. then assume worst case scenario * R1 is OK for inverse <=> CONDR1 .LT. REAL(N) * more conservative <=> CONDR1 .LT. SQRT(REAL(N)) @@ -1763,7 +1763,7 @@ ELSE * * .. ill-conditioned case: second QRF with pivoting -* Note that windowed pivoting would be equaly good +* Note that windowed pivoting would be equally good * numerically, and more run-time efficient. So, in * an optimal implementation, the next call to CGEQP3 * should be replaced with eg. CALL CGEQPX (ACM TOMS #782) @@ -1821,7 +1821,7 @@ * IF ( CONDR2 .GE. COND_OK ) THEN * .. save the Householder vectors used for Q3 -* (this overwrittes the copy of R2, as it will not be +* (this overwrites the copy of R2, as it will not be * needed in this branch, but it does not overwritte the * Huseholder vectors of Q2.). CALL CLACPY( 'U', NR, NR, V, LDV, CWORK(2*N+1), N ) @@ -2077,7 +2077,7 @@ * * This branch deploys a preconditioned Jacobi SVD with explicitly * accumulated rotations. It is included as optional, mainly for -* experimental purposes. It does perfom well, and can also be used. +* experimental purposes. It does perform well, and can also be used. * In this implementation, this branch will be automatically activated * if the condition number sigma_max(A) / sigma_min(A) is predicted * to be greater than the overflow threshold. This is because the diff --git a/lapack-netlib/SRC/cgelq.f b/lapack-netlib/SRC/cgelq.f index 909162ebc..c3b2238bf 100644 --- a/lapack-netlib/SRC/cgelq.f +++ b/lapack-netlib/SRC/cgelq.f @@ -1,3 +1,4 @@ +*> \brief \b CGELQ * * Definition: * =========== @@ -17,7 +18,17 @@ * ============= *> *> \verbatim -*> CGELQ computes a LQ factorization of an M-by-N matrix A. +*> +*> CGELQ computes an LQ factorization of a complex M-by-N matrix A: +*> +*> A = ( L 0 ) * Q +*> +*> where: +*> +*> Q is a N-by-N orthogonal matrix; +*> L is an lower-triangular M-by-M matrix; +*> 0 is a M-by-(N-M) zero matrix, if M < N. +*> *> \endverbatim * * Arguments: @@ -138,7 +149,7 @@ *> \verbatim *> *> These details are particular for this LAPACK implementation. Users should not -*> take them for granted. These details may change in the future, and are unlikely not +*> take them for granted. These details may change in the future, and are not likely *> true for another LAPACK implementation. These details are relevant if one wants *> to try to understand the code. They are not part of the interface. *> @@ -159,10 +170,10 @@ SUBROUTINE CGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK, $ INFO ) * -* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK computational routine (version 3.9.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. -- -* December 2016 +* November 2019 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N, TSIZE, LWORK diff --git a/lapack-netlib/SRC/cgelq2.f b/lapack-netlib/SRC/cgelq2.f index 9742d359b..3fab2c396 100644 --- a/lapack-netlib/SRC/cgelq2.f +++ b/lapack-netlib/SRC/cgelq2.f @@ -33,8 +33,16 @@ *> *> \verbatim *> -*> CGELQ2 computes an LQ factorization of a complex m by n matrix A: -*> A = L * Q. +*> CGELQ2 computes an LQ factorization of a complex m-by-n matrix A: +*> +*> A = ( L 0 ) * Q +*> +*> where: +*> +*> Q is a n-by-n orthogonal matrix; +*> L is an lower-triangular m-by-m matrix; +*> 0 is a m-by-(n-m) zero matrix, if m < n. +*> *> \endverbatim * * Arguments: @@ -96,7 +104,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date December 2016 +*> \date November 2019 * *> \ingroup complexGEcomputational * @@ -121,10 +129,10 @@ * ===================================================================== SUBROUTINE CGELQ2( M, N, A, LDA, TAU, WORK, INFO ) * -* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK computational routine (version 3.9.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* December 2016 +* November 2019 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N diff --git a/lapack-netlib/SRC/cgelqf.f b/lapack-netlib/SRC/cgelqf.f index 216630e88..030ac0b4d 100644 --- a/lapack-netlib/SRC/cgelqf.f +++ b/lapack-netlib/SRC/cgelqf.f @@ -34,7 +34,15 @@ *> \verbatim *> *> CGELQF computes an LQ factorization of a complex M-by-N matrix A: -*> A = L * Q. +*> +*> A = ( L 0 ) * Q +*> +*> where: +*> +*> Q is a N-by-N orthogonal matrix; +*> L is an lower-triangular M-by-M matrix; +*> 0 is a M-by-(N-M) zero matrix, if M < N. +*> *> \endverbatim * * Arguments: @@ -110,7 +118,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date December 2016 +*> \date November 2019 * *> \ingroup complexGEcomputational * @@ -135,10 +143,10 @@ * ===================================================================== SUBROUTINE CGELQF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) * -* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK computational routine (version 3.9.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* December 2016 +* November 2019 * * .. Scalar Arguments .. INTEGER INFO, LDA, LWORK, M, N diff --git a/lapack-netlib/SRC/cgelqt.f b/lapack-netlib/SRC/cgelqt.f index e151f10fe..f40db0b02 100644 --- a/lapack-netlib/SRC/cgelqt.f +++ b/lapack-netlib/SRC/cgelqt.f @@ -1,3 +1,4 @@ +*> \brief \b CGELQT * * Definition: * =========== diff --git a/lapack-netlib/SRC/cgelqt3.f b/lapack-netlib/SRC/cgelqt3.f index f64379722..80a9a9fc7 100644 --- a/lapack-netlib/SRC/cgelqt3.f +++ b/lapack-netlib/SRC/cgelqt3.f @@ -1,3 +1,5 @@ +*> \brief \b CGELQT3 +* * Definition: * =========== * diff --git a/lapack-netlib/SRC/cgemlq.f b/lapack-netlib/SRC/cgemlq.f index 2f44e7cfb..4e374077e 100644 --- a/lapack-netlib/SRC/cgemlq.f +++ b/lapack-netlib/SRC/cgemlq.f @@ -1,3 +1,4 @@ +*> \brief \b CGEMLQ * * Definition: * =========== @@ -143,7 +144,7 @@ *> \verbatim *> *> These details are particular for this LAPACK implementation. Users should not -*> take them for granted. These details may change in the future, and are unlikely not +*> take them for granted. These details may change in the future, and are not likely *> true for another LAPACK implementation. These details are relevant if one wants *> to try to understand the code. They are not part of the interface. *> diff --git a/lapack-netlib/SRC/cgemlqt.f b/lapack-netlib/SRC/cgemlqt.f index e35e421b1..66b186bff 100644 --- a/lapack-netlib/SRC/cgemlqt.f +++ b/lapack-netlib/SRC/cgemlqt.f @@ -1,3 +1,5 @@ +*> \brief \b CGEMLQT +* * Definition: * =========== * diff --git a/lapack-netlib/SRC/cgemqr.f b/lapack-netlib/SRC/cgemqr.f index a43d7be5b..54ab7aa74 100644 --- a/lapack-netlib/SRC/cgemqr.f +++ b/lapack-netlib/SRC/cgemqr.f @@ -1,3 +1,4 @@ +*> \brief \b CGEMQR * * Definition: * =========== @@ -144,7 +145,7 @@ *> \verbatim *> *> These details are particular for this LAPACK implementation. Users should not -*> take them for granted. These details may change in the future, and are unlikely not +*> take them for granted. These details may change in the future, and are not likely *> true for another LAPACK implementation. These details are relevant if one wants *> to try to understand the code. They are not part of the interface. *> diff --git a/lapack-netlib/SRC/cgeqr.f b/lapack-netlib/SRC/cgeqr.f index a00ef45c0..e0aea88b1 100644 --- a/lapack-netlib/SRC/cgeqr.f +++ b/lapack-netlib/SRC/cgeqr.f @@ -1,3 +1,4 @@ +*> \brief \b CGEQR * * Definition: * =========== @@ -17,7 +18,18 @@ * ============= *> *> \verbatim -*> CGEQR computes a QR factorization of an M-by-N matrix A. +*> +*> CGEQR computes a QR factorization of a complex M-by-N matrix A: +*> +*> A = Q * ( R ), +*> ( 0 ) +*> +*> where: +*> +*> Q is a M-by-M orthogonal matrix; +*> R is an upper-triangular N-by-N matrix; +*> 0 is a (M-N)-by-N zero matrix, if M > N. +*> *> \endverbatim * * Arguments: @@ -138,7 +150,7 @@ *> \verbatim *> *> These details are particular for this LAPACK implementation. Users should not -*> take them for granted. These details may change in the future, and are unlikely not +*> take them for granted. These details may change in the future, and are not likely *> true for another LAPACK implementation. These details are relevant if one wants *> to try to understand the code. They are not part of the interface. *> @@ -160,10 +172,10 @@ SUBROUTINE CGEQR( M, N, A, LDA, T, TSIZE, WORK, LWORK, $ INFO ) * -* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK computational routine (version 3.9.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. -- -* December 2016 +* November 2019 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N, TSIZE, LWORK diff --git a/lapack-netlib/SRC/cgeqr2.f b/lapack-netlib/SRC/cgeqr2.f index 1b2030b47..8cb2fa119 100644 --- a/lapack-netlib/SRC/cgeqr2.f +++ b/lapack-netlib/SRC/cgeqr2.f @@ -33,8 +33,17 @@ *> *> \verbatim *> -*> CGEQR2 computes a QR factorization of a complex m by n matrix A: -*> A = Q * R. +*> CGEQR2 computes a QR factorization of a complex m-by-n matrix A: +*> +*> A = Q * ( R ), +*> ( 0 ) +*> +*> where: +*> +*> Q is a m-by-m orthogonal matrix; +*> R is an upper-triangular n-by-n matrix; +*> 0 is a (m-n)-by-n zero matrix, if m > n. +*> *> \endverbatim * * Arguments: @@ -96,7 +105,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date December 2016 +*> \date November 2019 * *> \ingroup complexGEcomputational * @@ -121,10 +130,10 @@ * ===================================================================== SUBROUTINE CGEQR2( M, N, A, LDA, TAU, WORK, INFO ) * -* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK computational routine (version 3.9.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* December 2016 +* November 2019 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N diff --git a/lapack-netlib/SRC/cgeqr2p.f b/lapack-netlib/SRC/cgeqr2p.f index 3c64255d9..1e7b980df 100644 --- a/lapack-netlib/SRC/cgeqr2p.f +++ b/lapack-netlib/SRC/cgeqr2p.f @@ -33,8 +33,18 @@ *> *> \verbatim *> -*> CGEQR2P computes a QR factorization of a complex m by n matrix A: -*> A = Q * R. The diagonal entries of R are real and nonnegative. +*> CGEQR2P computes a QR factorization of a complex m-by-n matrix A: +*> +*> A = Q * ( R ), +*> ( 0 ) +*> +*> where: +*> +*> Q is a m-by-m orthogonal matrix; +*> R is an upper-triangular n-by-n matrix with nonnegative diagonal +*> entries; +*> 0 is a (m-n)-by-n zero matrix, if m > n. +*> *> \endverbatim * * Arguments: @@ -97,7 +107,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date December 2016 +*> \date November 2019 * *> \ingroup complexGEcomputational * @@ -124,10 +134,10 @@ * ===================================================================== SUBROUTINE CGEQR2P( M, N, A, LDA, TAU, WORK, INFO ) * -* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK computational routine (version 3.9.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* December 2016 +* November 2019 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N diff --git a/lapack-netlib/SRC/cgeqrf.f b/lapack-netlib/SRC/cgeqrf.f index 833384707..ff0c53f2f 100644 --- a/lapack-netlib/SRC/cgeqrf.f +++ b/lapack-netlib/SRC/cgeqrf.f @@ -34,7 +34,16 @@ *> \verbatim *> *> CGEQRF computes a QR factorization of a complex M-by-N matrix A: -*> A = Q * R. +*> +*> A = Q * ( R ), +*> ( 0 ) +*> +*> where: +*> +*> Q is a M-by-M orthogonal matrix; +*> R is an upper-triangular N-by-N matrix; +*> 0 is a (M-N)-by-N zero matrix, if M > N. +*> *> \endverbatim * * Arguments: @@ -111,7 +120,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date December 2016 +*> \date November 2019 * *> \ingroup complexGEcomputational * @@ -136,10 +145,10 @@ * ===================================================================== SUBROUTINE CGEQRF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) * -* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK computational routine (version 3.9.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* December 2016 +* November 2019 * * .. Scalar Arguments .. INTEGER INFO, LDA, LWORK, M, N diff --git a/lapack-netlib/SRC/cgeqrfp.f b/lapack-netlib/SRC/cgeqrfp.f index a56508b4e..9c29ac90b 100644 --- a/lapack-netlib/SRC/cgeqrfp.f +++ b/lapack-netlib/SRC/cgeqrfp.f @@ -33,8 +33,18 @@ *> *> \verbatim *> -*> CGEQRFP computes a QR factorization of a complex M-by-N matrix A: -*> A = Q * R. The diagonal entries of R are real and nonnegative. +*> CGEQR2P computes a QR factorization of a complex M-by-N matrix A: +*> +*> A = Q * ( R ), +*> ( 0 ) +*> +*> where: +*> +*> Q is a M-by-M orthogonal matrix; +*> R is an upper-triangular N-by-N matrix with nonnegative diagonal +*> entries; +*> 0 is a (M-N)-by-N zero matrix, if M > N. +*> *> \endverbatim * * Arguments: @@ -112,7 +122,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date December 2016 +*> \date November 2019 * *> \ingroup complexGEcomputational * @@ -139,10 +149,10 @@ * ===================================================================== SUBROUTINE CGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO ) * -* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK computational routine (version 3.9.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* December 2016 +* November 2019 * * .. Scalar Arguments .. INTEGER INFO, LDA, LWORK, M, N diff --git a/lapack-netlib/SRC/cgerfsx.f b/lapack-netlib/SRC/cgerfsx.f index 7b72f9c9a..a6e24ae4f 100644 --- a/lapack-netlib/SRC/cgerfsx.f +++ b/lapack-netlib/SRC/cgerfsx.f @@ -74,7 +74,7 @@ *> Specifies the form of the system of equations: *> = 'N': A * X = B (No transpose) *> = 'T': A**T * X = B (Transpose) -*> = 'C': A**H * X = B (Conjugate transpose = Transpose) +*> = 'C': A**H * X = B (Conjugate transpose) *> \endverbatim *> *> \param[in] EQUED @@ -283,7 +283,7 @@ *> information as described below. There currently are up to three *> pieces of information returned for each right-hand side. If *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then -*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most +*> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most *> the first (:,N_ERR_BNDS) entries are returned. *> *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith @@ -319,14 +319,14 @@ *> \param[in] NPARAMS *> \verbatim *> NPARAMS is INTEGER -*> Specifies the number of parameters set in PARAMS. If .LE. 0, the +*> Specifies the number of parameters set in PARAMS. If <= 0, the *> PARAMS array is never referenced and default values are used. *> \endverbatim *> *> \param[in,out] PARAMS *> \verbatim *> PARAMS is REAL array, dimension NPARAMS -*> Specifies algorithm parameters. If an entry is .LT. 0.0, then +*> Specifies algorithm parameters. If an entry is < 0.0, then *> that entry will be filled with default value used for that *> parameter. Only positions up to NPARAMS are accessed; defaults *> are used for higher-numbered parameters. @@ -334,9 +334,9 @@ *> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative *> refinement or not. *> Default: 1.0 -*> = 0.0 : No refinement is performed, and no error bounds are +*> = 0.0: No refinement is performed, and no error bounds are *> computed. -*> = 1.0 : Use the double-precision refinement algorithm, +*> = 1.0: Use the double-precision refinement algorithm, *> possibly with doubled-single computations if the *> compilation environment does not support DOUBLE *> PRECISION. diff --git a/lapack-netlib/SRC/cgesc2.f b/lapack-netlib/SRC/cgesc2.f index c0b91107e..6f45a09a6 100644 --- a/lapack-netlib/SRC/cgesc2.f +++ b/lapack-netlib/SRC/cgesc2.f @@ -91,7 +91,7 @@ *> \verbatim *> SCALE is REAL *> On exit, SCALE contains the scale factor. SCALE is chosen -*> 0 <= SCALE <= 1 to prevent owerflow in the solution. +*> 0 <= SCALE <= 1 to prevent overflow in the solution. *> \endverbatim * * Authors: diff --git a/lapack-netlib/SRC/cgesvdq.f b/lapack-netlib/SRC/cgesvdq.f new file mode 100644 index 000000000..77c883dde --- /dev/null +++ b/lapack-netlib/SRC/cgesvdq.f @@ -0,0 +1,1391 @@ +*> \brief CGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CGESVDQ + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA, +* S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK, +* CWORK, LCWORK, RWORK, LRWORK, INFO ) +* +* .. Scalar Arguments .. +* IMPLICIT NONE +* CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV +* INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LCWORK, LRWORK, +* INFO +* .. +* .. Array Arguments .. +* COMPLEX A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( * ) +* REAL S( * ), RWORK( * ) +* INTEGER IWORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CGESVDQ computes the singular value decomposition (SVD) of a complex +*> M-by-N matrix A, where M >= N. The SVD of A is written as +*> [++] [xx] [x0] [xx] +*> A = U * SIGMA * V^*, [++] = [xx] * [ox] * [xx] +*> [++] [xx] +*> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal +*> matrix, and V is an N-by-N unitary matrix. The diagonal elements +*> of SIGMA are the singular values of A. The columns of U and V are the +*> left and the right singular vectors of A, respectively. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] JOBA +*> \verbatim +*> JOBA is CHARACTER*1 +*> Specifies the level of accuracy in the computed SVD +*> = 'A' The requested accuracy corresponds to having the backward +*> error bounded by || delta A ||_F <= f(m,n) * EPS * || A ||_F, +*> where EPS = SLAMCH('Epsilon'). This authorises CGESVDQ to +*> truncate the computed triangular factor in a rank revealing +*> QR factorization whenever the truncated part is below the +*> threshold of the order of EPS * ||A||_F. This is aggressive +*> truncation level. +*> = 'M' Similarly as with 'A', but the truncation is more gentle: it +*> is allowed only when there is a drop on the diagonal of the +*> triangular factor in the QR factorization. This is medium +*> truncation level. +*> = 'H' High accuracy requested. No numerical rank determination based +*> on the rank revealing QR factorization is attempted. +*> = 'E' Same as 'H', and in addition the condition number of column +*> scaled A is estimated and returned in RWORK(1). +*> N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1) +*> \endverbatim +*> +*> \param[in] JOBP +*> \verbatim +*> JOBP is CHARACTER*1 +*> = 'P' The rows of A are ordered in decreasing order with respect to +*> ||A(i,:)||_\infty. This enhances numerical accuracy at the cost +*> of extra data movement. Recommended for numerical robustness. +*> = 'N' No row pivoting. +*> \endverbatim +*> +*> \param[in] JOBR +*> \verbatim +*> JOBR is CHARACTER*1 +*> = 'T' After the initial pivoted QR factorization, CGESVD is applied to +*> the adjoint R**H of the computed triangular factor R. This involves +*> some extra data movement (matrix transpositions). Useful for +*> experiments, research and development. +*> = 'N' The triangular factor R is given as input to CGESVD. This may be +*> preferred as it involves less data movement. +*> \endverbatim +*> +*> \param[in] JOBU +*> \verbatim +*> JOBU is CHARACTER*1 +*> = 'A' All M left singular vectors are computed and returned in the +*> matrix U. See the description of U. +*> = 'S' or 'U' N = min(M,N) left singular vectors are computed and returned +*> in the matrix U. See the description of U. +*> = 'R' Numerical rank NUMRANK is determined and only NUMRANK left singular +*> vectors are computed and returned in the matrix U. +*> = 'F' The N left singular vectors are returned in factored form as the +*> product of the Q factor from the initial QR factorization and the +*> N left singular vectors of (R**H , 0)**H. If row pivoting is used, +*> then the necessary information on the row pivoting is stored in +*> IWORK(N+1:N+M-1). +*> = 'N' The left singular vectors are not computed. +*> \endverbatim +*> +*> \param[in] JOBV +*> \verbatim +*> JOBV is CHARACTER*1 +*> = 'A', 'V' All N right singular vectors are computed and returned in +*> the matrix V. +*> = 'R' Numerical rank NUMRANK is determined and only NUMRANK right singular +*> vectors are computed and returned in the matrix V. This option is +*> allowed only if JOBU = 'R' or JOBU = 'N'; otherwise it is illegal. +*> = 'N' The right singular vectors are not computed. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the input matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the input matrix A. M >= N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX array of dimensions LDA x N +*> On entry, the input matrix A. +*> On exit, if JOBU .NE. 'N' or JOBV .NE. 'N', the lower triangle of A contains +*> the Householder vectors as stored by CGEQP3. If JOBU = 'F', these Householder +*> vectors together with CWORK(1:N) can be used to restore the Q factors from +*> the initial pivoted QR factorization of A. See the description of U. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER. +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[out] S +*> \verbatim +*> S is REAL array of dimension N. +*> The singular values of A, ordered so that S(i) >= S(i+1). +*> \endverbatim +*> +*> \param[out] U +*> \verbatim +*> U is COMPLEX array, dimension +*> LDU x M if JOBU = 'A'; see the description of LDU. In this case, +*> on exit, U contains the M left singular vectors. +*> LDU x N if JOBU = 'S', 'U', 'R' ; see the description of LDU. In this +*> case, U contains the leading N or the leading NUMRANK left singular vectors. +*> LDU x N if JOBU = 'F' ; see the description of LDU. In this case U +*> contains N x N unitary matrix that can be used to form the left +*> singular vectors. +*> If JOBU = 'N', U is not referenced. +*> \endverbatim +*> +*> \param[in] LDU +*> \verbatim +*> LDU is INTEGER. +*> The leading dimension of the array U. +*> If JOBU = 'A', 'S', 'U', 'R', LDU >= max(1,M). +*> If JOBU = 'F', LDU >= max(1,N). +*> Otherwise, LDU >= 1. +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is COMPLEX array, dimension +*> LDV x N if JOBV = 'A', 'V', 'R' or if JOBA = 'E' . +*> If JOBV = 'A', or 'V', V contains the N-by-N unitary matrix V**H; +*> If JOBV = 'R', V contains the first NUMRANK rows of V**H (the right +*> singular vectors, stored rowwise, of the NUMRANK largest singular values). +*> If JOBV = 'N' and JOBA = 'E', V is used as a workspace. +*> If JOBV = 'N', and JOBA.NE.'E', V is not referenced. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. +*> If JOBV = 'A', 'V', 'R', or JOBA = 'E', LDV >= max(1,N). +*> Otherwise, LDV >= 1. +*> \endverbatim +*> +*> \param[out] NUMRANK +*> \verbatim +*> NUMRANK is INTEGER +*> NUMRANK is the numerical rank first determined after the rank +*> revealing QR factorization, following the strategy specified by the +*> value of JOBA. If JOBV = 'R' and JOBU = 'R', only NUMRANK +*> leading singular values and vectors are then requested in the call +*> of CGESVD. The final value of NUMRANK might be further reduced if +*> some singular values are computed as zeros. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (max(1, LIWORK)). +*> On exit, IWORK(1:N) contains column pivoting permutation of the +*> rank revealing QR factorization. +*> If JOBP = 'P', IWORK(N+1:N+M-1) contains the indices of the sequence +*> of row swaps used in row pivoting. These can be used to restore the +*> left singular vectors in the case JOBU = 'F'. +*> +*> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0, +*> LIWORK(1) returns the minimal LIWORK. +*> \endverbatim +*> +*> \param[in] LIWORK +*> \verbatim +*> LIWORK is INTEGER +*> The dimension of the array IWORK. +*> LIWORK >= N + M - 1, if JOBP = 'P'; +*> LIWORK >= N if JOBP = 'N'. +*> +*> If LIWORK = -1, then a workspace query is assumed; the routine +*> only calculates and returns the optimal and minimal sizes +*> for the CWORK, IWORK, and RWORK arrays, and no error +*> message related to LCWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] CWORK +*> \verbatim +*> CWORK is COMPLEX array, dimension (max(2, LCWORK)), used as a workspace. +*> On exit, if, on entry, LCWORK.NE.-1, CWORK(1:N) contains parameters +*> needed to recover the Q factor from the QR factorization computed by +*> CGEQP3. +*> +*> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0, +*> CWORK(1) returns the optimal LCWORK, and +*> CWORK(2) returns the minimal LCWORK. +*> \endverbatim +*> +*> \param[in,out] LCWORK +*> \verbatim +*> LCWORK is INTEGER +*> The dimension of the array CWORK. It is determined as follows: +*> Let LWQP3 = N+1, LWCON = 2*N, and let +*> LWUNQ = { MAX( N, 1 ), if JOBU = 'R', 'S', or 'U' +*> { MAX( M, 1 ), if JOBU = 'A' +*> LWSVD = MAX( 3*N, 1 ) +*> LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 3*(N/2), 1 ), LWUNLQ = MAX( N, 1 ), +*> LWQRF = MAX( N/2, 1 ), LWUNQ2 = MAX( N, 1 ) +*> Then the minimal value of LCWORK is: +*> = MAX( N + LWQP3, LWSVD ) if only the singular values are needed; +*> = MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values are needed, +*> and a scaled condition estimate requested; +*> +*> = N + MAX( LWQP3, LWSVD, LWUNQ ) if the singular values and the left +*> singular vectors are requested; +*> = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the singular values and the left +*> singular vectors are requested, and also +*> a scaled condition estimate requested; +*> +*> = N + MAX( LWQP3, LWSVD ) if the singular values and the right +*> singular vectors are requested; +*> = N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and the right +*> singular vectors are requested, and also +*> a scaled condition etimate requested; +*> +*> = N + MAX( LWQP3, LWSVD, LWUNQ ) if the full SVD is requested with JOBV = 'R'; +*> independent of JOBR; +*> = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the full SVD is requested, +*> JOBV = 'R' and, also a scaled condition +*> estimate requested; independent of JOBR; +*> = MAX( N + MAX( LWQP3, LWSVD, LWUNQ ), +*> N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ) ) if the +*> full SVD is requested with JOBV = 'A' or 'V', and +*> JOBR ='N' +*> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ), +*> N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ ) ) +*> if the full SVD is requested with JOBV = 'A' or 'V', and +*> JOBR ='N', and also a scaled condition number estimate +*> requested. +*> = MAX( N + MAX( LWQP3, LWSVD, LWUNQ ), +*> N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) ) if the +*> full SVD is requested with JOBV = 'A', 'V', and JOBR ='T' +*> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ), +*> N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) ) +*> if the full SVD is requested with JOBV = 'A', 'V' and +*> JOBR ='T', and also a scaled condition number estimate +*> requested. +*> Finally, LCWORK must be at least two: LCWORK = MAX( 2, LCWORK ). +*> +*> If LCWORK = -1, then a workspace query is assumed; the routine +*> only calculates and returns the optimal and minimal sizes +*> for the CWORK, IWORK, and RWORK arrays, and no error +*> message related to LCWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (max(1, LRWORK)). +*> On exit, +*> 1. If JOBA = 'E', RWORK(1) contains an estimate of the condition +*> number of column scaled A. If A = C * D where D is diagonal and C +*> has unit columns in the Euclidean norm, then, assuming full column rank, +*> N^(-1/4) * RWORK(1) <= ||pinv(C)||_2 <= N^(1/4) * RWORK(1). +*> Otherwise, RWORK(1) = -1. +*> 2. RWORK(2) contains the number of singular values computed as +*> exact zeros in CGESVD applied to the upper triangular or trapeziodal +*> R (from the initial QR factorization). In case of early exit (no call to +*> CGESVD, such as in the case of zero matrix) RWORK(2) = -1. +*> +*> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0, +*> RWORK(1) returns the minimal LRWORK. +*> \endverbatim +*> +*> \param[in] LRWORK +*> \verbatim +*> LRWORK is INTEGER. +*> The dimension of the array RWORK. +*> If JOBP ='P', then LRWORK >= MAX(2, M, 5*N); +*> Otherwise, LRWORK >= MAX(2, 5*N). +*> +*> If LRWORK = -1, then a workspace query is assumed; the routine +*> only calculates and returns the optimal and minimal sizes +*> for the CWORK, IWORK, and RWORK arrays, and no error +*> message related to LCWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit. +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> > 0: if CBDSQR did not converge, INFO specifies how many superdiagonals +*> of an intermediate bidiagonal form B (computed in CGESVD) did not +*> converge to zero. +*> \endverbatim +* +*> \par Further Details: +* ======================== +*> +*> \verbatim +*> +*> 1. The data movement (matrix transpose) is coded using simple nested +*> DO-loops because BLAS and LAPACK do not provide corresponding subroutines. +*> Those DO-loops are easily identified in this source code - by the CONTINUE +*> statements labeled with 11**. In an optimized version of this code, the +*> nested DO loops should be replaced with calls to an optimized subroutine. +*> 2. This code scales A by 1/SQRT(M) if the largest ABS(A(i,j)) could cause +*> column norm overflow. This is the minial precaution and it is left to the +*> SVD routine (CGESVD) to do its own preemptive scaling if potential over- +*> or underflows are detected. To avoid repeated scanning of the array A, +*> an optimal implementation would do all necessary scaling before calling +*> CGESVD and the scaling in CGESVD can be switched off. +*> 3. Other comments related to code optimization are given in comments in the +*> code, enlosed in [[double brackets]]. +*> \endverbatim +* +*> \par Bugs, examples and comments +* =========================== +* +*> \verbatim +*> Please report all bugs and send interesting examples and/or comments to +*> drmac@math.hr. Thank you. +*> \endverbatim +* +*> \par References +* =============== +* +*> \verbatim +*> [1] Zlatko Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for +*> Computing the SVD with High Accuracy. ACM Trans. Math. Softw. +*> 44(1): 11:1-11:30 (2017) +*> +*> SIGMA library, xGESVDQ section updated February 2016. +*> Developed and coded by Zlatko Drmac, Department of Mathematics +*> University of Zagreb, Croatia, drmac@math.hr +*> \endverbatim +* +* +*> \par Contributors: +* ================== +*> +*> \verbatim +*> Developed and coded by Zlatko Drmac, Department of Mathematics +*> University of Zagreb, Croatia, drmac@math.hr +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2018 +* +*> \ingroup complexGEsing +* +* ===================================================================== + SUBROUTINE CGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA, + $ S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK, + $ CWORK, LCWORK, RWORK, LRWORK, INFO ) +* .. Scalar Arguments .. + IMPLICIT NONE + CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV + INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LCWORK, LRWORK, + $ INFO +* .. +* .. Array Arguments .. + COMPLEX A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( * ) + REAL S( * ), RWORK( * ) + INTEGER IWORK( * ) +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) + COMPLEX CZERO, CONE + PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ), CONE = ( 1.0E0, 0.0E0 ) ) +* .. +* .. Local Scalars .. + INTEGER IERR, NR, N1, OPTRATIO, p, q + INTEGER LWCON, LWQP3, LWRK_CGELQF, LWRK_CGESVD, LWRK_CGESVD2, + $ LWRK_CGEQP3, LWRK_CGEQRF, LWRK_CUNMLQ, LWRK_CUNMQR, + $ LWRK_CUNMQR2, LWLQF, LWQRF, LWSVD, LWSVD2, LWUNQ, + $ LWUNQ2, LWUNLQ, MINWRK, MINWRK2, OPTWRK, OPTWRK2, + $ IMINWRK, RMINWRK + LOGICAL ACCLA, ACCLM, ACCLH, ASCALED, CONDA, DNTWU, DNTWV, + $ LQUERY, LSVC0, LSVEC, ROWPRM, RSVEC, RTRANS, WNTUA, + $ WNTUF, WNTUR, WNTUS, WNTVA, WNTVR + REAL BIG, EPSLN, RTMP, SCONDA, SFMIN + COMPLEX CTMP +* .. +* .. Local Arrays + COMPLEX CDUMMY(1) + REAL RDUMMY(1) +* .. +* .. External Subroutines (BLAS, LAPACK) + EXTERNAL CGELQF, CGEQP3, CGEQRF, CGESVD, CLACPY, CLAPMT, + $ CLASCL, CLASET, CLASWP, CSSCAL, SLASET, SLASCL, + $ CPOCON, CUNMLQ, CUNMQR, XERBLA +* .. +* .. External Functions (BLAS, LAPACK) + LOGICAL LSAME + INTEGER ISAMAX + REAL CLANGE, SCNRM2, SLAMCH + EXTERNAL CLANGE, LSAME, ISAMAX, SCNRM2, SLAMCH +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, CONJG, MAX, MIN, REAL, SQRT +* .. +* .. Executable Statements .. +* +* Test the input arguments +* + WNTUS = LSAME( JOBU, 'S' ) .OR. LSAME( JOBU, 'U' ) + WNTUR = LSAME( JOBU, 'R' ) + WNTUA = LSAME( JOBU, 'A' ) + WNTUF = LSAME( JOBU, 'F' ) + LSVC0 = WNTUS .OR. WNTUR .OR. WNTUA + LSVEC = LSVC0 .OR. WNTUF + DNTWU = LSAME( JOBU, 'N' ) +* + WNTVR = LSAME( JOBV, 'R' ) + WNTVA = LSAME( JOBV, 'A' ) .OR. LSAME( JOBV, 'V' ) + RSVEC = WNTVR .OR. WNTVA + DNTWV = LSAME( JOBV, 'N' ) +* + ACCLA = LSAME( JOBA, 'A' ) + ACCLM = LSAME( JOBA, 'M' ) + CONDA = LSAME( JOBA, 'E' ) + ACCLH = LSAME( JOBA, 'H' ) .OR. CONDA +* + ROWPRM = LSAME( JOBP, 'P' ) + RTRANS = LSAME( JOBR, 'T' ) +* + IF ( ROWPRM ) THEN + IMINWRK = MAX( 1, N + M - 1 ) + RMINWRK = MAX( 2, M, 5*N ) + ELSE + IMINWRK = MAX( 1, N ) + RMINWRK = MAX( 2, 5*N ) + END IF + LQUERY = (LIWORK .EQ. -1 .OR. LCWORK .EQ. -1 .OR. LRWORK .EQ. -1) + INFO = 0 + IF ( .NOT. ( ACCLA .OR. ACCLM .OR. ACCLH ) ) THEN + INFO = -1 + ELSE IF ( .NOT.( ROWPRM .OR. LSAME( JOBP, 'N' ) ) ) THEN + INFO = -2 + ELSE IF ( .NOT.( RTRANS .OR. LSAME( JOBR, 'N' ) ) ) THEN + INFO = -3 + ELSE IF ( .NOT.( LSVEC .OR. DNTWU ) ) THEN + INFO = -4 + ELSE IF ( WNTUR .AND. WNTVA ) THEN + INFO = -5 + ELSE IF ( .NOT.( RSVEC .OR. DNTWV )) THEN + INFO = -5 + ELSE IF ( M.LT.0 ) THEN + INFO = -6 + ELSE IF ( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN + INFO = -7 + ELSE IF ( LDA.LT.MAX( 1, M ) ) THEN + INFO = -9 + ELSE IF ( LDU.LT.1 .OR. ( LSVC0 .AND. LDU.LT.M ) .OR. + $ ( WNTUF .AND. LDU.LT.N ) ) THEN + INFO = -12 + ELSE IF ( LDV.LT.1 .OR. ( RSVEC .AND. LDV.LT.N ) .OR. + $ ( CONDA .AND. LDV.LT.N ) ) THEN + INFO = -14 + ELSE IF ( LIWORK .LT. IMINWRK .AND. .NOT. LQUERY ) THEN + INFO = -17 + END IF +* +* + IF ( INFO .EQ. 0 ) THEN +* +* Compute workspace +* .. compute the minimal and the optimal workspace lengths +* [[The expressions for computing the minimal and the optimal +* values of LCWORK are written with a lot of redundancy and +* can be simplified. However, this detailed form is easier for +* maintenance and modifications of the code.]] +* +* .. minimal workspace length for CGEQP3 of an M x N matrix + LWQP3 = N+1 +* .. minimal workspace length for CUNMQR to build left singular vectors + IF ( WNTUS .OR. WNTUR ) THEN + LWUNQ = MAX( N , 1 ) + ELSE IF ( WNTUA ) THEN + LWUNQ = MAX( M , 1 ) + END IF +* .. minimal workspace length for CPOCON of an N x N matrix + LWCON = 2 * N +* .. CGESVD of an N x N matrix + LWSVD = MAX( 3 * N, 1 ) + IF ( LQUERY ) THEN + CALL CGEQP3( M, N, A, LDA, IWORK, CDUMMY, CDUMMY, -1, + $ RDUMMY, IERR ) + LWRK_CGEQP3 = INT( CDUMMY(1) ) + IF ( WNTUS .OR. WNTUR ) THEN + CALL CUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U, + $ LDU, CDUMMY, -1, IERR ) + LWRK_CUNMQR = INT( CDUMMY(1) ) + ELSE IF ( WNTUA ) THEN + CALL CUNMQR( 'L', 'N', M, M, N, A, LDA, CDUMMY, U, + $ LDU, CDUMMY, -1, IERR ) + LWRK_CUNMQR = INT( CDUMMY(1) ) + ELSE + LWRK_CUNMQR = 0 + END IF + END IF + MINWRK = 2 + OPTWRK = 2 + IF ( .NOT. (LSVEC .OR. RSVEC )) THEN +* .. minimal and optimal sizes of the complex workspace if +* only the singular values are requested + IF ( CONDA ) THEN + MINWRK = MAX( N+LWQP3, LWCON, LWSVD ) + ELSE + MINWRK = MAX( N+LWQP3, LWSVD ) + END IF + IF ( LQUERY ) THEN + CALL CGESVD( 'N', 'N', N, N, A, LDA, S, U, LDU, + $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) + LWRK_CGESVD = INT( CDUMMY(1) ) + IF ( CONDA ) THEN + OPTWRK = MAX( N+LWRK_CGEQP3, N+LWCON, LWRK_CGESVD ) + ELSE + OPTWRK = MAX( N+LWRK_CGEQP3, LWRK_CGESVD ) + END IF + END IF + ELSE IF ( LSVEC .AND. (.NOT.RSVEC) ) THEN +* .. minimal and optimal sizes of the complex workspace if the +* singular values and the left singular vectors are requested + IF ( CONDA ) THEN + MINWRK = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) + ELSE + MINWRK = N + MAX( LWQP3, LWSVD, LWUNQ ) + END IF + IF ( LQUERY ) THEN + IF ( RTRANS ) THEN + CALL CGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU, + $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) + ELSE + CALL CGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU, + $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) + END IF + LWRK_CGESVD = INT( CDUMMY(1) ) + IF ( CONDA ) THEN + OPTWRK = N + MAX( LWRK_CGEQP3, LWCON, LWRK_CGESVD, + $ LWRK_CUNMQR ) + ELSE + OPTWRK = N + MAX( LWRK_CGEQP3, LWRK_CGESVD, + $ LWRK_CUNMQR ) + END IF + END IF + ELSE IF ( RSVEC .AND. (.NOT.LSVEC) ) THEN +* .. minimal and optimal sizes of the complex workspace if the +* singular values and the right singular vectors are requested + IF ( CONDA ) THEN + MINWRK = N + MAX( LWQP3, LWCON, LWSVD ) + ELSE + MINWRK = N + MAX( LWQP3, LWSVD ) + END IF + IF ( LQUERY ) THEN + IF ( RTRANS ) THEN + CALL CGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU, + $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) + ELSE + CALL CGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU, + $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) + END IF + LWRK_CGESVD = INT( CDUMMY(1) ) + IF ( CONDA ) THEN + OPTWRK = N + MAX( LWRK_CGEQP3, LWCON, LWRK_CGESVD ) + ELSE + OPTWRK = N + MAX( LWRK_CGEQP3, LWRK_CGESVD ) + END IF + END IF + ELSE +* .. minimal and optimal sizes of the complex workspace if the +* full SVD is requested + IF ( RTRANS ) THEN + MINWRK = MAX( LWQP3, LWSVD, LWUNQ ) + IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON ) + MINWRK = MINWRK + N + IF ( WNTVA ) THEN +* .. minimal workspace length for N x N/2 CGEQRF + LWQRF = MAX( N/2, 1 ) +* .. minimal workspace lengt for N/2 x N/2 CGESVD + LWSVD2 = MAX( 3 * (N/2), 1 ) + LWUNQ2 = MAX( N, 1 ) + MINWRK2 = MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, + $ N/2+LWUNQ2, LWUNQ ) + IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON ) + MINWRK2 = N + MINWRK2 + MINWRK = MAX( MINWRK, MINWRK2 ) + END IF + ELSE + MINWRK = MAX( LWQP3, LWSVD, LWUNQ ) + IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON ) + MINWRK = MINWRK + N + IF ( WNTVA ) THEN +* .. minimal workspace length for N/2 x N CGELQF + LWLQF = MAX( N/2, 1 ) + LWSVD2 = MAX( 3 * (N/2), 1 ) + LWUNLQ = MAX( N , 1 ) + MINWRK2 = MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, + $ N/2+LWUNLQ, LWUNQ ) + IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON ) + MINWRK2 = N + MINWRK2 + MINWRK = MAX( MINWRK, MINWRK2 ) + END IF + END IF + IF ( LQUERY ) THEN + IF ( RTRANS ) THEN + CALL CGESVD( 'O', 'A', N, N, A, LDA, S, U, LDU, + $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) + LWRK_CGESVD = INT( CDUMMY(1) ) + OPTWRK = MAX(LWRK_CGEQP3,LWRK_CGESVD,LWRK_CUNMQR) + IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON ) + OPTWRK = N + OPTWRK + IF ( WNTVA ) THEN + CALL CGEQRF(N,N/2,U,LDU,CDUMMY,CDUMMY,-1,IERR) + LWRK_CGEQRF = INT( CDUMMY(1) ) + CALL CGESVD( 'S', 'O', N/2,N/2, V,LDV, S, U,LDU, + $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) + LWRK_CGESVD2 = INT( CDUMMY(1) ) + CALL CUNMQR( 'R', 'C', N, N, N/2, U, LDU, CDUMMY, + $ V, LDV, CDUMMY, -1, IERR ) + LWRK_CUNMQR2 = INT( CDUMMY(1) ) + OPTWRK2 = MAX( LWRK_CGEQP3, N/2+LWRK_CGEQRF, + $ N/2+LWRK_CGESVD2, N/2+LWRK_CUNMQR2 ) + IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON ) + OPTWRK2 = N + OPTWRK2 + OPTWRK = MAX( OPTWRK, OPTWRK2 ) + END IF + ELSE + CALL CGESVD( 'S', 'O', N, N, A, LDA, S, U, LDU, + $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) + LWRK_CGESVD = INT( CDUMMY(1) ) + OPTWRK = MAX(LWRK_CGEQP3,LWRK_CGESVD,LWRK_CUNMQR) + IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON ) + OPTWRK = N + OPTWRK + IF ( WNTVA ) THEN + CALL CGELQF(N/2,N,U,LDU,CDUMMY,CDUMMY,-1,IERR) + LWRK_CGELQF = INT( CDUMMY(1) ) + CALL CGESVD( 'S','O', N/2,N/2, V, LDV, S, U, LDU, + $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) + LWRK_CGESVD2 = INT( CDUMMY(1) ) + CALL CUNMLQ( 'R', 'N', N, N, N/2, U, LDU, CDUMMY, + $ V, LDV, CDUMMY,-1,IERR ) + LWRK_CUNMLQ = INT( CDUMMY(1) ) + OPTWRK2 = MAX( LWRK_CGEQP3, N/2+LWRK_CGELQF, + $ N/2+LWRK_CGESVD2, N/2+LWRK_CUNMLQ ) + IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON ) + OPTWRK2 = N + OPTWRK2 + OPTWRK = MAX( OPTWRK, OPTWRK2 ) + END IF + END IF + END IF + END IF +* + MINWRK = MAX( 2, MINWRK ) + OPTWRK = MAX( 2, OPTWRK ) + IF ( LCWORK .LT. MINWRK .AND. (.NOT.LQUERY) ) INFO = -19 +* + END IF +* + IF (INFO .EQ. 0 .AND. LRWORK .LT. RMINWRK .AND. .NOT. LQUERY) THEN + INFO = -21 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CGESVDQ', -INFO ) + RETURN + ELSE IF ( LQUERY ) THEN +* +* Return optimal workspace +* + IWORK(1) = IMINWRK + CWORK(1) = OPTWRK + CWORK(2) = MINWRK + RWORK(1) = RMINWRK + RETURN + END IF +* +* Quick return if the matrix is void. +* + IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) ) THEN +* .. all output is void. + RETURN + END IF +* + BIG = SLAMCH('O') + ASCALED = .FALSE. + IF ( ROWPRM ) THEN +* .. reordering the rows in decreasing sequence in the +* ell-infinity norm - this enhances numerical robustness in +* the case of differently scaled rows. + DO 1904 p = 1, M +* RWORK(p) = ABS( A(p,ICAMAX(N,A(p,1),LDA)) ) +* [[CLANGE will return NaN if an entry of the p-th row is Nan]] + RWORK(p) = CLANGE( 'M', 1, N, A(p,1), LDA, RDUMMY ) +* .. check for NaN's and Inf's + IF ( ( RWORK(p) .NE. RWORK(p) ) .OR. + $ ( (RWORK(p)*ZERO) .NE. ZERO ) ) THEN + INFO = - 8 + CALL XERBLA( 'CGESVDQ', -INFO ) + RETURN + END IF + 1904 CONTINUE + DO 1952 p = 1, M - 1 + q = ISAMAX( M-p+1, RWORK(p), 1 ) + p - 1 + IWORK(N+p) = q + IF ( p .NE. q ) THEN + RTMP = RWORK(p) + RWORK(p) = RWORK(q) + RWORK(q) = RTMP + END IF + 1952 CONTINUE +* + IF ( RWORK(1) .EQ. ZERO ) THEN +* Quick return: A is the M x N zero matrix. + NUMRANK = 0 + CALL SLASET( 'G', N, 1, ZERO, ZERO, S, N ) + IF ( WNTUS ) CALL CLASET('G', M, N, CZERO, CONE, U, LDU) + IF ( WNTUA ) CALL CLASET('G', M, M, CZERO, CONE, U, LDU) + IF ( WNTVA ) CALL CLASET('G', N, N, CZERO, CONE, V, LDV) + IF ( WNTUF ) THEN + CALL CLASET( 'G', N, 1, CZERO, CZERO, CWORK, N ) + CALL CLASET( 'G', M, N, CZERO, CONE, U, LDU ) + END IF + DO 5001 p = 1, N + IWORK(p) = p + 5001 CONTINUE + IF ( ROWPRM ) THEN + DO 5002 p = N + 1, N + M - 1 + IWORK(p) = p - N + 5002 CONTINUE + END IF + IF ( CONDA ) RWORK(1) = -1 + RWORK(2) = -1 + RETURN + END IF +* + IF ( RWORK(1) .GT. BIG / SQRT(REAL(M)) ) THEN +* .. to prevent overflow in the QR factorization, scale the +* matrix by 1/sqrt(M) if too large entry detected + CALL CLASCL('G',0,0,SQRT(REAL(M)),ONE, M,N, A,LDA, IERR) + ASCALED = .TRUE. + END IF + CALL CLASWP( N, A, LDA, 1, M-1, IWORK(N+1), 1 ) + END IF +* +* .. At this stage, preemptive scaling is done only to avoid column +* norms overflows during the QR factorization. The SVD procedure should +* have its own scaling to save the singular values from overflows and +* underflows. That depends on the SVD procedure. +* + IF ( .NOT.ROWPRM ) THEN + RTMP = CLANGE( 'M', M, N, A, LDA, RWORK ) + IF ( ( RTMP .NE. RTMP ) .OR. + $ ( (RTMP*ZERO) .NE. ZERO ) ) THEN + INFO = - 8 + CALL XERBLA( 'CGESVDQ', -INFO ) + RETURN + END IF + IF ( RTMP .GT. BIG / SQRT(REAL(M)) ) THEN +* .. to prevent overflow in the QR factorization, scale the +* matrix by 1/sqrt(M) if too large entry detected + CALL CLASCL('G',0,0, SQRT(REAL(M)),ONE, M,N, A,LDA, IERR) + ASCALED = .TRUE. + END IF + END IF +* +* .. QR factorization with column pivoting +* +* A * P = Q * [ R ] +* [ 0 ] +* + DO 1963 p = 1, N +* .. all columns are free columns + IWORK(p) = 0 + 1963 CONTINUE + CALL CGEQP3( M, N, A, LDA, IWORK, CWORK, CWORK(N+1), LCWORK-N, + $ RWORK, IERR ) +* +* If the user requested accuracy level allows truncation in the +* computed upper triangular factor, the matrix R is examined and, +* if possible, replaced with its leading upper trapezoidal part. +* + EPSLN = SLAMCH('E') + SFMIN = SLAMCH('S') +* SMALL = SFMIN / EPSLN + NR = N +* + IF ( ACCLA ) THEN +* +* Standard absolute error bound suffices. All sigma_i with +* sigma_i < N*EPS*||A||_F are flushed to zero. This is an +* aggressive enforcement of lower numerical rank by introducing a +* backward error of the order of N*EPS*||A||_F. + NR = 1 + RTMP = SQRT(REAL(N))*EPSLN + DO 3001 p = 2, N + IF ( ABS(A(p,p)) .LT. (RTMP*ABS(A(1,1))) ) GO TO 3002 + NR = NR + 1 + 3001 CONTINUE + 3002 CONTINUE +* + ELSEIF ( ACCLM ) THEN +* .. similarly as above, only slightly more gentle (less aggressive). +* Sudden drop on the diagonal of R is used as the criterion for being +* close-to-rank-deficient. The threshold is set to EPSLN=SLAMCH('E'). +* [[This can be made more flexible by replacing this hard-coded value +* with a user specified threshold.]] Also, the values that underflow +* will be truncated. + NR = 1 + DO 3401 p = 2, N + IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR. + $ ( ABS(A(p,p)) .LT. SFMIN ) ) GO TO 3402 + NR = NR + 1 + 3401 CONTINUE + 3402 CONTINUE +* + ELSE +* .. RRQR not authorized to determine numerical rank except in the +* obvious case of zero pivots. +* .. inspect R for exact zeros on the diagonal; +* R(i,i)=0 => R(i:N,i:N)=0. + NR = 1 + DO 3501 p = 2, N + IF ( ABS(A(p,p)) .EQ. ZERO ) GO TO 3502 + NR = NR + 1 + 3501 CONTINUE + 3502 CONTINUE +* + IF ( CONDA ) THEN +* Estimate the scaled condition number of A. Use the fact that it is +* the same as the scaled condition number of R. +* .. V is used as workspace + CALL CLACPY( 'U', N, N, A, LDA, V, LDV ) +* Only the leading NR x NR submatrix of the triangular factor +* is considered. Only if NR=N will this give a reliable error +* bound. However, even for NR < N, this can be used on an +* expert level and obtain useful information in the sense of +* perturbation theory. + DO 3053 p = 1, NR + RTMP = SCNRM2( p, V(1,p), 1 ) + CALL CSSCAL( p, ONE/RTMP, V(1,p), 1 ) + 3053 CONTINUE + IF ( .NOT. ( LSVEC .OR. RSVEC ) ) THEN + CALL CPOCON( 'U', NR, V, LDV, ONE, RTMP, + $ CWORK, RWORK, IERR ) + ELSE + CALL CPOCON( 'U', NR, V, LDV, ONE, RTMP, + $ CWORK(N+1), RWORK, IERR ) + END IF + SCONDA = ONE / SQRT(RTMP) +* For NR=N, SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1), +* N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA +* See the reference [1] for more details. + END IF +* + ENDIF +* + IF ( WNTUR ) THEN + N1 = NR + ELSE IF ( WNTUS .OR. WNTUF) THEN + N1 = N + ELSE IF ( WNTUA ) THEN + N1 = M + END IF +* + IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN +*....................................................................... +* .. only the singular values are requested +*....................................................................... + IF ( RTRANS ) THEN +* +* .. compute the singular values of R**H = [A](1:NR,1:N)**H +* .. set the lower triangle of [A] to [A](1:NR,1:N)**H and +* the upper triangle of [A] to zero. + DO 1146 p = 1, MIN( N, NR ) + A(p,p) = CONJG(A(p,p)) + DO 1147 q = p + 1, N + A(q,p) = CONJG(A(p,q)) + IF ( q .LE. NR ) A(p,q) = CZERO + 1147 CONTINUE + 1146 CONTINUE +* + CALL CGESVD( 'N', 'N', N, NR, A, LDA, S, U, LDU, + $ V, LDV, CWORK, LCWORK, RWORK, INFO ) +* + ELSE +* +* .. compute the singular values of R = [A](1:NR,1:N) +* + IF ( NR .GT. 1 ) + $ CALL CLASET( 'L', NR-1,NR-1, CZERO,CZERO, A(2,1), LDA ) + CALL CGESVD( 'N', 'N', NR, N, A, LDA, S, U, LDU, + $ V, LDV, CWORK, LCWORK, RWORK, INFO ) +* + END IF +* + ELSE IF ( LSVEC .AND. ( .NOT. RSVEC) ) THEN +*....................................................................... +* .. the singular values and the left singular vectors requested +*......................................................................."""""""" + IF ( RTRANS ) THEN +* .. apply CGESVD to R**H +* .. copy R**H into [U] and overwrite [U] with the right singular +* vectors of R + DO 1192 p = 1, NR + DO 1193 q = p, N + U(q,p) = CONJG(A(p,q)) + 1193 CONTINUE + 1192 CONTINUE + IF ( NR .GT. 1 ) + $ CALL CLASET( 'U', NR-1,NR-1, CZERO,CZERO, U(1,2), LDU ) +* .. the left singular vectors not computed, the NR right singular +* vectors overwrite [U](1:NR,1:NR) as conjugate transposed. These +* will be pre-multiplied by Q to build the left singular vectors of A. + CALL CGESVD( 'N', 'O', N, NR, U, LDU, S, U, LDU, + $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO ) +* + DO 1119 p = 1, NR + U(p,p) = CONJG(U(p,p)) + DO 1120 q = p + 1, NR + CTMP = CONJG(U(q,p)) + U(q,p) = CONJG(U(p,q)) + U(p,q) = CTMP + 1120 CONTINUE + 1119 CONTINUE +* + ELSE +* .. apply CGESVD to R +* .. copy R into [U] and overwrite [U] with the left singular vectors + CALL CLACPY( 'U', NR, N, A, LDA, U, LDU ) + IF ( NR .GT. 1 ) + $ CALL CLASET( 'L', NR-1, NR-1, CZERO, CZERO, U(2,1), LDU ) +* .. the right singular vectors not computed, the NR left singular +* vectors overwrite [U](1:NR,1:NR) + CALL CGESVD( 'O', 'N', NR, N, U, LDU, S, U, LDU, + $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO ) +* .. now [U](1:NR,1:NR) contains the NR left singular vectors of +* R. These will be pre-multiplied by Q to build the left singular +* vectors of A. + END IF +* +* .. assemble the left singular vector matrix U of dimensions +* (M x NR) or (M x N) or (M x M). + IF ( ( NR .LT. M ) .AND. ( .NOT.WNTUF ) ) THEN + CALL CLASET('A', M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU) + IF ( NR .LT. N1 ) THEN + CALL CLASET( 'A',NR,N1-NR,CZERO,CZERO,U(1,NR+1), LDU ) + CALL CLASET( 'A',M-NR,N1-NR,CZERO,CONE, + $ U(NR+1,NR+1), LDU ) + END IF + END IF +* +* The Q matrix from the first QRF is built into the left singular +* vectors matrix U. +* + IF ( .NOT.WNTUF ) + $ CALL CUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U, + $ LDU, CWORK(N+1), LCWORK-N, IERR ) + IF ( ROWPRM .AND. .NOT.WNTUF ) + $ CALL CLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 ) +* + ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN +*....................................................................... +* .. the singular values and the right singular vectors requested +*....................................................................... + IF ( RTRANS ) THEN +* .. apply CGESVD to R**H +* .. copy R**H into V and overwrite V with the left singular vectors + DO 1165 p = 1, NR + DO 1166 q = p, N + V(q,p) = CONJG(A(p,q)) + 1166 CONTINUE + 1165 CONTINUE + IF ( NR .GT. 1 ) + $ CALL CLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV ) +* .. the left singular vectors of R**H overwrite V, the right singular +* vectors not computed + IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN + CALL CGESVD( 'O', 'N', N, NR, V, LDV, S, U, LDU, + $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO ) +* + DO 1121 p = 1, NR + V(p,p) = CONJG(V(p,p)) + DO 1122 q = p + 1, NR + CTMP = CONJG(V(q,p)) + V(q,p) = CONJG(V(p,q)) + V(p,q) = CTMP + 1122 CONTINUE + 1121 CONTINUE +* + IF ( NR .LT. N ) THEN + DO 1103 p = 1, NR + DO 1104 q = NR + 1, N + V(p,q) = CONJG(V(q,p)) + 1104 CONTINUE + 1103 CONTINUE + END IF + CALL CLAPMT( .FALSE., NR, N, V, LDV, IWORK ) + ELSE +* .. need all N right singular vectors and NR < N +* [!] This is simple implementation that augments [V](1:N,1:NR) +* by padding a zero block. In the case NR << N, a more efficient +* way is to first use the QR factorization. For more details +* how to implement this, see the " FULL SVD " branch. + CALL CLASET('G', N, N-NR, CZERO, CZERO, V(1,NR+1), LDV) + CALL CGESVD( 'O', 'N', N, N, V, LDV, S, U, LDU, + $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO ) +* + DO 1123 p = 1, N + V(p,p) = CONJG(V(p,p)) + DO 1124 q = p + 1, N + CTMP = CONJG(V(q,p)) + V(q,p) = CONJG(V(p,q)) + V(p,q) = CTMP + 1124 CONTINUE + 1123 CONTINUE + CALL CLAPMT( .FALSE., N, N, V, LDV, IWORK ) + END IF +* + ELSE +* .. aply CGESVD to R +* .. copy R into V and overwrite V with the right singular vectors + CALL CLACPY( 'U', NR, N, A, LDA, V, LDV ) + IF ( NR .GT. 1 ) + $ CALL CLASET( 'L', NR-1, NR-1, CZERO, CZERO, V(2,1), LDV ) +* .. the right singular vectors overwrite V, the NR left singular +* vectors stored in U(1:NR,1:NR) + IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN + CALL CGESVD( 'N', 'O', NR, N, V, LDV, S, U, LDU, + $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO ) + CALL CLAPMT( .FALSE., NR, N, V, LDV, IWORK ) +* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H + ELSE +* .. need all N right singular vectors and NR < N +* [!] This is simple implementation that augments [V](1:NR,1:N) +* by padding a zero block. In the case NR << N, a more efficient +* way is to first use the LQ factorization. For more details +* how to implement this, see the " FULL SVD " branch. + CALL CLASET('G', N-NR, N, CZERO,CZERO, V(NR+1,1), LDV) + CALL CGESVD( 'N', 'O', N, N, V, LDV, S, U, LDU, + $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO ) + CALL CLAPMT( .FALSE., N, N, V, LDV, IWORK ) + END IF +* .. now [V] contains the adjoint of the matrix of the right singular +* vectors of A. + END IF +* + ELSE +*....................................................................... +* .. FULL SVD requested +*....................................................................... + IF ( RTRANS ) THEN +* +* .. apply CGESVD to R**H [[this option is left for R&D&T]] +* + IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN +* .. copy R**H into [V] and overwrite [V] with the left singular +* vectors of R**H + DO 1168 p = 1, NR + DO 1169 q = p, N + V(q,p) = CONJG(A(p,q)) + 1169 CONTINUE + 1168 CONTINUE + IF ( NR .GT. 1 ) + $ CALL CLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV ) +* +* .. the left singular vectors of R**H overwrite [V], the NR right +* singular vectors of R**H stored in [U](1:NR,1:NR) as conjugate +* transposed + CALL CGESVD( 'O', 'A', N, NR, V, LDV, S, V, LDV, + $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO ) +* .. assemble V + DO 1115 p = 1, NR + V(p,p) = CONJG(V(p,p)) + DO 1116 q = p + 1, NR + CTMP = CONJG(V(q,p)) + V(q,p) = CONJG(V(p,q)) + V(p,q) = CTMP + 1116 CONTINUE + 1115 CONTINUE + IF ( NR .LT. N ) THEN + DO 1101 p = 1, NR + DO 1102 q = NR+1, N + V(p,q) = CONJG(V(q,p)) + 1102 CONTINUE + 1101 CONTINUE + END IF + CALL CLAPMT( .FALSE., NR, N, V, LDV, IWORK ) +* + DO 1117 p = 1, NR + U(p,p) = CONJG(U(p,p)) + DO 1118 q = p + 1, NR + CTMP = CONJG(U(q,p)) + U(q,p) = CONJG(U(p,q)) + U(p,q) = CTMP + 1118 CONTINUE + 1117 CONTINUE +* + IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN + CALL CLASET('A', M-NR,NR, CZERO,CZERO, U(NR+1,1), LDU) + IF ( NR .LT. N1 ) THEN + CALL CLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU) + CALL CLASET( 'A',M-NR,N1-NR,CZERO,CONE, + $ U(NR+1,NR+1), LDU ) + END IF + END IF +* + ELSE +* .. need all N right singular vectors and NR < N +* .. copy R**H into [V] and overwrite [V] with the left singular +* vectors of R**H +* [[The optimal ratio N/NR for using QRF instead of padding +* with zeros. Here hard coded to 2; it must be at least +* two due to work space constraints.]] +* OPTRATIO = ILAENV(6, 'CGESVD', 'S' // 'O', NR,N,0,0) +* OPTRATIO = MAX( OPTRATIO, 2 ) + OPTRATIO = 2 + IF ( OPTRATIO*NR .GT. N ) THEN + DO 1198 p = 1, NR + DO 1199 q = p, N + V(q,p) = CONJG(A(p,q)) + 1199 CONTINUE + 1198 CONTINUE + IF ( NR .GT. 1 ) + $ CALL CLASET('U',NR-1,NR-1, CZERO,CZERO, V(1,2),LDV) +* + CALL CLASET('A',N,N-NR,CZERO,CZERO,V(1,NR+1),LDV) + CALL CGESVD( 'O', 'A', N, N, V, LDV, S, V, LDV, + $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO ) +* + DO 1113 p = 1, N + V(p,p) = CONJG(V(p,p)) + DO 1114 q = p + 1, N + CTMP = CONJG(V(q,p)) + V(q,p) = CONJG(V(p,q)) + V(p,q) = CTMP + 1114 CONTINUE + 1113 CONTINUE + CALL CLAPMT( .FALSE., N, N, V, LDV, IWORK ) +* .. assemble the left singular vector matrix U of dimensions +* (M x N1), i.e. (M x N) or (M x M). +* + DO 1111 p = 1, N + U(p,p) = CONJG(U(p,p)) + DO 1112 q = p + 1, N + CTMP = CONJG(U(q,p)) + U(q,p) = CONJG(U(p,q)) + U(p,q) = CTMP + 1112 CONTINUE + 1111 CONTINUE +* + IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) THEN + CALL CLASET('A',M-N,N,CZERO,CZERO,U(N+1,1),LDU) + IF ( N .LT. N1 ) THEN + CALL CLASET('A',N,N1-N,CZERO,CZERO,U(1,N+1),LDU) + CALL CLASET('A',M-N,N1-N,CZERO,CONE, + $ U(N+1,N+1), LDU ) + END IF + END IF + ELSE +* .. copy R**H into [U] and overwrite [U] with the right +* singular vectors of R + DO 1196 p = 1, NR + DO 1197 q = p, N + U(q,NR+p) = CONJG(A(p,q)) + 1197 CONTINUE + 1196 CONTINUE + IF ( NR .GT. 1 ) + $ CALL CLASET('U',NR-1,NR-1,CZERO,CZERO,U(1,NR+2),LDU) + CALL CGEQRF( N, NR, U(1,NR+1), LDU, CWORK(N+1), + $ CWORK(N+NR+1), LCWORK-N-NR, IERR ) + DO 1143 p = 1, NR + DO 1144 q = 1, N + V(q,p) = CONJG(U(p,NR+q)) + 1144 CONTINUE + 1143 CONTINUE + CALL CLASET('U',NR-1,NR-1,CZERO,CZERO,V(1,2),LDV) + CALL CGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU, + $ V,LDV, CWORK(N+NR+1),LCWORK-N-NR,RWORK, INFO ) + CALL CLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV) + CALL CLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV) + CALL CLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV) + CALL CUNMQR('R','C', N, N, NR, U(1,NR+1), LDU, + $ CWORK(N+1),V,LDV,CWORK(N+NR+1),LCWORK-N-NR,IERR) + CALL CLAPMT( .FALSE., N, N, V, LDV, IWORK ) +* .. assemble the left singular vector matrix U of dimensions +* (M x NR) or (M x N) or (M x M). + IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN + CALL CLASET('A',M-NR,NR,CZERO,CZERO,U(NR+1,1),LDU) + IF ( NR .LT. N1 ) THEN + CALL CLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU) + CALL CLASET( 'A',M-NR,N1-NR,CZERO,CONE, + $ U(NR+1,NR+1),LDU) + END IF + END IF + END IF + END IF +* + ELSE +* +* .. apply CGESVD to R [[this is the recommended option]] +* + IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN +* .. copy R into [V] and overwrite V with the right singular vectors + CALL CLACPY( 'U', NR, N, A, LDA, V, LDV ) + IF ( NR .GT. 1 ) + $ CALL CLASET( 'L', NR-1,NR-1, CZERO,CZERO, V(2,1), LDV ) +* .. the right singular vectors of R overwrite [V], the NR left +* singular vectors of R stored in [U](1:NR,1:NR) + CALL CGESVD( 'S', 'O', NR, N, V, LDV, S, U, LDU, + $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO ) + CALL CLAPMT( .FALSE., NR, N, V, LDV, IWORK ) +* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H +* .. assemble the left singular vector matrix U of dimensions +* (M x NR) or (M x N) or (M x M). + IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN + CALL CLASET('A', M-NR,NR, CZERO,CZERO, U(NR+1,1), LDU) + IF ( NR .LT. N1 ) THEN + CALL CLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU) + CALL CLASET( 'A',M-NR,N1-NR,CZERO,CONE, + $ U(NR+1,NR+1), LDU ) + END IF + END IF +* + ELSE +* .. need all N right singular vectors and NR < N +* .. the requested number of the left singular vectors +* is then N1 (N or M) +* [[The optimal ratio N/NR for using LQ instead of padding +* with zeros. Here hard coded to 2; it must be at least +* two due to work space constraints.]] +* OPTRATIO = ILAENV(6, 'CGESVD', 'S' // 'O', NR,N,0,0) +* OPTRATIO = MAX( OPTRATIO, 2 ) + OPTRATIO = 2 + IF ( OPTRATIO * NR .GT. N ) THEN + CALL CLACPY( 'U', NR, N, A, LDA, V, LDV ) + IF ( NR .GT. 1 ) + $ CALL CLASET('L', NR-1,NR-1, CZERO,CZERO, V(2,1),LDV) +* .. the right singular vectors of R overwrite [V], the NR left +* singular vectors of R stored in [U](1:NR,1:NR) + CALL CLASET('A', N-NR,N, CZERO,CZERO, V(NR+1,1),LDV) + CALL CGESVD( 'S', 'O', N, N, V, LDV, S, U, LDU, + $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO ) + CALL CLAPMT( .FALSE., N, N, V, LDV, IWORK ) +* .. now [V] contains the adjoint of the matrix of the right +* singular vectors of A. The leading N left singular vectors +* are in [U](1:N,1:N) +* .. assemble the left singular vector matrix U of dimensions +* (M x N1), i.e. (M x N) or (M x M). + IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) THEN + CALL CLASET('A',M-N,N,CZERO,CZERO,U(N+1,1),LDU) + IF ( N .LT. N1 ) THEN + CALL CLASET('A',N,N1-N,CZERO,CZERO,U(1,N+1),LDU) + CALL CLASET( 'A',M-N,N1-N,CZERO,CONE, + $ U(N+1,N+1), LDU ) + END IF + END IF + ELSE + CALL CLACPY( 'U', NR, N, A, LDA, U(NR+1,1), LDU ) + IF ( NR .GT. 1 ) + $ CALL CLASET('L',NR-1,NR-1,CZERO,CZERO,U(NR+2,1),LDU) + CALL CGELQF( NR, N, U(NR+1,1), LDU, CWORK(N+1), + $ CWORK(N+NR+1), LCWORK-N-NR, IERR ) + CALL CLACPY('L',NR,NR,U(NR+1,1),LDU,V,LDV) + IF ( NR .GT. 1 ) + $ CALL CLASET('U',NR-1,NR-1,CZERO,CZERO,V(1,2),LDV) + CALL CGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU, + $ V, LDV, CWORK(N+NR+1), LCWORK-N-NR, RWORK, INFO ) + CALL CLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV) + CALL CLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV) + CALL CLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV) + CALL CUNMLQ('R','N',N,N,NR,U(NR+1,1),LDU,CWORK(N+1), + $ V, LDV, CWORK(N+NR+1),LCWORK-N-NR,IERR) + CALL CLAPMT( .FALSE., N, N, V, LDV, IWORK ) +* .. assemble the left singular vector matrix U of dimensions +* (M x NR) or (M x N) or (M x M). + IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN + CALL CLASET('A',M-NR,NR,CZERO,CZERO,U(NR+1,1),LDU) + IF ( NR .LT. N1 ) THEN + CALL CLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU) + CALL CLASET( 'A',M-NR,N1-NR,CZERO,CONE, + $ U(NR+1,NR+1), LDU ) + END IF + END IF + END IF + END IF +* .. end of the "R**H or R" branch + END IF +* +* The Q matrix from the first QRF is built into the left singular +* vectors matrix U. +* + IF ( .NOT. WNTUF ) + $ CALL CUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U, + $ LDU, CWORK(N+1), LCWORK-N, IERR ) + IF ( ROWPRM .AND. .NOT.WNTUF ) + $ CALL CLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 ) +* +* ... end of the "full SVD" branch + END IF +* +* Check whether some singular values are returned as zeros, e.g. +* due to underflow, and update the numerical rank. + p = NR + DO 4001 q = p, 1, -1 + IF ( S(q) .GT. ZERO ) GO TO 4002 + NR = NR - 1 + 4001 CONTINUE + 4002 CONTINUE +* +* .. if numerical rank deficiency is detected, the truncated +* singular values are set to zero. + IF ( NR .LT. N ) CALL SLASET( 'G', N-NR,1, ZERO,ZERO, S(NR+1), N ) +* .. undo scaling; this may cause overflow in the largest singular +* values. + IF ( ASCALED ) + $ CALL SLASCL( 'G',0,0, ONE,SQRT(REAL(M)), NR,1, S, N, IERR ) + IF ( CONDA ) RWORK(1) = SCONDA + RWORK(2) = p - NR +* .. p-NR is the number of singular values that are computed as +* exact zeros in CGESVD() applied to the (possibly truncated) +* full row rank triangular (trapezoidal) factor of A. + NUMRANK = NR +* + RETURN +* +* End of CGESVDQ +* + END