lapackblas.h File Reference
#include <cmath>
Include dependency graph for lapackblas.h:
This graph shows which files directly or indirectly include this file:
Go to the source code of this file.
Defines | |
| #define | TRUE_ (1) |
| #define | FALSE_ (0) |
| #define | abs(x) ((x) >= 0 ? (x) : -(x)) |
| #define | dabs(x) (doublereal)abs(x) |
| #define | f2cmin(a, b) ((a) <= (b) ? (a) : (b)) |
| #define | f2cmax(a, b) ((a) >= (b) ? (a) : (b)) |
| #define | df2cmin(a, b) (doublereal)f2cmin(a,b) |
| #define | df2cmax(a, b) (doublereal)f2cmax(a,b) |
Typedefs | |
| typedef int | integer |
| typedef float | real |
| typedef double | doublereal |
| typedef int | logical |
| typedef short | flag |
| typedef short | ftnlen |
| typedef short | ftnint |
Functions | |
| int | s_cat (char *, const char **, integer *, integer *, ftnlen) |
| integer | s_cmp (char *, const char *, ftnlen, ftnlen) |
| void | s_copy (char *a, const char *b, ftnlen la, ftnlen lb) |
| double | pow_ri (real *ap, integer *bp) |
| double | r_sign (real *a, real *b) |
| integer | ieeeck_ (integer *ispec, real *zero, real *one) |
| integer | ilaenv_ (integer *ispec, const char *name__, const char *opts, integer *n1, integer *n2, integer *n3, integer *n4, ftnlen name_len, ftnlen opts_len) |
| real | drand (void) |
| logical | lsame_ (const char *ca, const char *cb) |
| int | saxpy_ (integer *n, real *sa, real *sx, integer *incx, real *sy, integer *incy) |
| int | scopy_ (integer *n, real *sx, integer *incx, real *sy, integer *incy) |
| doublereal | sdot_ (integer *n, real *sx, integer *incx, real *sy, integer *incy) |
| int | sgemm_ (const char *transa, const char *transb, integer *m, integer *n, integer *k, real *alpha, real *a, integer *lda, real *b, integer *ldb, real *beta, real *c__, integer *ldc) |
| int | sgemv_ (const char *trans, integer *m, integer *n, real *alpha, real *a, integer *lda, real *x, integer *incx, real *beta, real *y, integer *incy) |
| int | sger_ (integer *m, integer *n, real *alpha, real *x, integer *incx, real *y, integer *incy, real *a, integer *lda) |
| int | slae2_ (real *a, real *b, real *c__, real *rt1, real *rt2) |
| int | slaev2_ (real *a, real *b, real *c__, real *rt1, real *rt2, real *cs1, real *sn1) |
| doublereal | slamch_ (const char *cmach) |
| doublereal | slanst_ (const char *norm, integer *n, real *d__, real *e) |
| doublereal | slansy_ (const char *norm, char *uplo, integer *n, real *a, integer *lda, real *work) |
| doublereal | slapy2_ (real *x, real *y) |
| int | slarfb_ (const char *side, const char *trans, const char *direct, const char *storev, integer *m, integer *n, integer *k, real *v, integer *ldv, real *t, integer *ldt, real *c__, integer *ldc, real *work, integer *ldwork) |
| int | slarf_ (const char *side, integer *m, integer *n, real *v, integer *incv, real *tau, real *c__, integer *ldc, real *work) |
| int | slarfg_ (integer *n, real *alpha, real *x, integer *incx, real *tau) |
| int | slarft_ (const char *direct, const char *storev, integer *n, integer *k, real *v, integer *ldv, real *tau, real *t, integer *ldt) |
| int | slartg_ (real *f, real *g, real *cs, real *sn, real *r__) |
| int | slascl_ (const char *type__, integer *kl, integer *ku, real *cfrom, real *cto, integer *m, integer *n, real *a, integer *lda, integer *info) |
| int | slaset_ (const char *uplo, integer *m, integer *n, real *alpha, real *beta, real *a, integer *lda) |
| int | slasr_ (const char *side, const char *pivot, const char *direct, integer *m, integer *n, real *c__, real *s, real *a, integer *lda) |
| int | slasrt_ (const char *id, integer *n, real *d__, integer *info) |
| int | slassq_ (integer *n, real *x, integer *incx, real *scale, real *sumsq) |
| int | slatrd_ (char *uplo, integer *n, integer *nb, real *a, integer *lda, real *e, real *tau, real *w, integer *ldw) |
| doublereal | snrm2_ (integer *n, real *x, integer *incx) |
| int | srot_ (integer *n, real *sx, integer *incx, real *sy, integer *incy, real *c__, real *s) |
| int | sorg2l_ (integer *m, integer *n, integer *k, real *a, integer *lda, real *tau, real *work, integer *info) |
| int | sorg2r_ (integer *m, integer *n, integer *k, real *a, integer *lda, real *tau, real *work, integer *info) |
| int | sorgql_ (integer *m, integer *n, integer *k, real *a, integer *lda, real *tau, real *work, integer *lwork, integer *info) |
| int | sorgqr_ (integer *m, integer *n, integer *k, real *a, integer *lda, real *tau, real *work, integer *lwork, integer *info) |
| int | sorgtr_ (char *uplo, integer *n, real *a, integer *lda, real *tau, real *work, integer *lwork, integer *info) |
| int | sscal_ (integer *n, real *sa, real *sx, integer *incx) |
| int | ssteqr_ (const char *compz, integer *n, real *d__, real *e, real *z__, integer *ldz, real *work, integer *info) |
| int | ssterf_ (integer *n, real *d__, real *e, integer *info) |
| int | sswap_ (integer *n, real *sx, integer *incx, real *sy, integer *incy) |
| int | ssyev_ (char *jobz, char *uplo, integer *n, real *a, integer *lda, real *w, real *work, integer *lwork, integer *info) |
| int | ssymv_ (const char *uplo, integer *n, real *alpha, real *a, integer *lda, real *x, integer *incx, real *beta, real *y, integer *incy) |
| int | ssyr2_ (char *uplo, integer *n, real *alpha, real *x, integer *incx, real *y, integer *incy, real *a, integer *lda) |
| int | ssyr2k_ (char *uplo, const char *trans, integer *n, integer *k, real *alpha, real *a, integer *lda, real *b, integer *ldb, real *beta, real *c__, integer *ldc) |
| int | ssytd2_ (char *uplo, integer *n, real *a, integer *lda, real *d__, real *e, real *tau, integer *info) |
| int | ssytrd_ (char *uplo, integer *n, real *a, integer *lda, real *d__, real *e, real *tau, real *work, integer *lwork, integer *info) |
| int | strmm_ (const char *side, const char *uplo, const char *transa, const char *diag, integer *m, integer *n, real *alpha, real *a, integer *lda, real *b, integer *ldb) |
| int | strmv_ (const char *uplo, const char *trans, const char *diag, integer *n, real *a, integer *lda, real *x, integer *incx) |
| int | xerbla_ (const char *srname, integer *info) |
| int | sstedc_ (const char *compz, integer *n, real *d__, real *e, real *z__, integer *ldz, real *work, integer *lwork, integer *iwork, integer *liwork, integer *info) |
| int | sstevd_ (char *jobz, integer *n, real *d__, real *e, real *z__, integer *ldz, real *work, integer *lwork, integer *iwork, integer *liwork, integer *info) |
| int | slaeda_ (integer *n, integer *tlvls, integer *curlvl, integer *curpbm, integer *prmptr, integer *perm, integer *givptr, integer *givcol, real *givnum, real *q, integer *qptr, real *z__, real *ztemp, integer *info) |
| int | slaed0_ (integer *icompq, integer *qsiz, integer *n, real *d__, real *e, real *q, integer *ldq, real *qstore, integer *ldqs, real *work, integer *iwork, integer *info) |
| int | slaed1_ (integer *n, real *d__, real *q, integer *ldq, integer *indxq, real *rho, integer *cutpnt, real *work, integer *iwork, integer *info) |
| int | slaed2_ (integer *k, integer *n, integer *n1, real *d__, real *q, integer *ldq, integer *indxq, real *rho, real *z__, real *dlamda, real *w, real *q2, integer *indx, integer *indxc, integer *indxp, integer *coltyp, integer *info) |
| int | slaed3_ (integer *k, integer *n, integer *n1, real *d__, real *q, integer *ldq, real *rho, real *dlamda, real *q2, integer *indx, integer *ctot, real *w, real *s, integer *info) |
| int | slaed4_ (integer *n, integer *i__, real *d__, real *z__, real *delta, real *rho, real *dlam, integer *info) |
| int | slaed5_ (integer *i__, real *d__, real *z__, real *delta, real *rho, real *dlam) |
| int | slaed6_ (integer *kniter, logical *orgati, real *rho, real *d__, real *z__, real *finit, real *tau, integer *info) |
| int | slaed7_ (integer *icompq, integer *n, integer *qsiz, integer *tlvls, integer *curlvl, integer *curpbm, real *d__, real *q, integer *ldq, integer *indxq, real *rho, integer *cutpnt, real *qstore, integer *qptr, integer *prmptr, integer *perm, integer *givptr, integer *givcol, real *givnum, real *work, integer *iwork, integer *info) |
| int | slaed8_ (integer *icompq, integer *k, integer *n, integer *qsiz, real *d__, real *q, integer *ldq, integer *indxq, real *rho, integer *cutpnt, real *z__, real *dlamda, real *q2, integer *ldq2, real *w, integer *perm, integer *givptr, integer *givcol, real *givnum, integer *indxp, integer *indx, integer *info) |
| int | slaed9_ (integer *k, integer *kstart, integer *kstop, integer *n, real *d__, real *q, integer *ldq, real *rho, real *dlamda, real *w, real *s, integer *lds, integer *info) |
| int | slacpy_ (const char *uplo, integer *m, integer *n, real *a, integer *lda, real *b, integer *ldb) |
| int | slamrg_ (integer *n1, integer *n2, real *a, integer *strd1, integer *strd2, integer *index) |
| int | sorgl2_ (integer *m, integer *n, integer *k, real *a, integer *lda, real *tau, real *work, integer *info) |
| int | sgesvd_ (char *jobu, char *jobvt, integer *m, integer *n, real *a, integer *lda, real *s, real *u, integer *ldu, real *vt, integer *ldvt, real *work, integer *lwork, integer *info) |
| int | sorglq_ (integer *m, integer *n, integer *k, real *a, integer *lda, real *tau, real *work, integer *lwork, integer *info) |
| doublereal | slange_ (const char *norm, integer *m, integer *n, real *a, integer *lda, real *work) |
| int | sgebrd_ (integer *m, integer *n, real *a, integer *lda, real *d__, real *e, real *tauq, real *taup, real *work, integer *lwork, integer *info) |
| int | sgebd2_ (integer *m, integer *n, real *a, integer *lda, real *d__, real *e, real *tauq, real *taup, real *work, integer *info) |
| int | sormbr_ (const char *vect, const char *side, const char *trans, integer *m, integer *n, integer *k, real *a, integer *lda, real *tau, real *c__, integer *ldc, real *work, integer *lwork, integer *info) |
| int | sgelqf_ (integer *m, integer *n, real *a, integer *lda, real *tau, real *work, integer *lwork, integer *info) |
| int | sormlq_ (const char *side, const char *trans, integer *m, integer *n, integer *k, real *a, integer *lda, real *tau, real *c__, integer *ldc, real *work, integer *lwork, integer *info) |
| int | sormqr_ (const char *side, const char *trans, integer *m, integer *n, integer *k, real *a, integer *lda, real *tau, real *c__, integer *ldc, real *work, integer *lwork, integer *info) |
| int | sgelq2_ (integer *m, integer *n, real *a, integer *lda, real *tau, real *work, integer *info) |
| int | sbdsqr_ (const char *uplo, integer *n, integer *ncvt, integer *nru, integer *ncc, real *d__, real *e, real *vt, integer *ldvt, real *u, integer *ldu, real *c__, integer *ldc, real *work, integer *info) |
| int | sgeqrf_ (integer *m, integer *n, real *a, integer *lda, real *tau, real *work, integer *lwork, integer *info) |
| int | sorml2_ (const char *side, const char *trans, integer *m, integer *n, integer *k, real *a, integer *lda, real *tau, real *c__, integer *ldc, real *work, integer *info) |
| int | slabrd_ (integer *m, integer *n, integer *nb, real *a, integer *lda, real *d__, real *e, real *tauq, real *taup, real *x, integer *ldx, real *y, integer *ldy) |
| int | sgeqr2_ (integer *m, integer *n, real *a, integer *lda, real *tau, real *work, integer *info) |
| int | sorm2r_ (const char *side, const char *trans, integer *m, integer *n, integer *k, real *a, integer *lda, real *tau, real *c__, integer *ldc, real *work, integer *info) |
| int | sorgbr_ (const char *vect, integer *m, integer *n, integer *k, real *a, integer *lda, real *tau, real *work, integer *lwork, integer *info) |
| int | slasq1_ (integer *n, real *d__, real *e, real *work, integer *info) |
| int | slasq2_ (integer *n, real *z__, integer *info) |
| int | slasq3_ (integer *i0, integer *n0, real *z__, integer *pp, real *dmin__, real *sigma, real *desig, real *qmax, integer *nfail, integer *iter, integer *ndiv, logical *ieee) |
| int | slasq4_ (integer *i0, integer *n0, real *z__, integer *pp, integer *n0in, real *dmin__, real *dmin1, real *dmin2, real *dn, real *dn1, real *dn2, real *tau, integer *ttype) |
| int | slasq5_ (integer *i0, integer *n0, real *z__, integer *pp, real *tau, real *dmin__, real *dmin1, real *dmin2, real *dn, real *dnm1, real *dnm2, logical *ieee) |
| int | slasq6_ (integer *i0, integer *n0, real *z__, integer *pp, real *dmin__, real *dmin1, real *dmin2, real *dn, real *dnm1, real *dnm2) |
| int | slasv2_ (real *f, real *g, real *h__, real *ssmin, real *ssmax, real *snr, real *csr, real *snl, real *csl) |
| int | slas2_ (real *f, real *g, real *h__, real *ssmin, real *ssmax) |
Define Documentation
| #define abs | ( | x | ) | ((x) >= 0 ? (x) : -(x)) |
Definition at line 46 of file lapackblas.h.
Referenced by EMAN::EMData::absi(), addnod_(), aprq2d(), EMAN::Util::cl1(), EMAN::Util::cml_line_insino(), EMAN::Util::cml_line_insino_all(), crlist_(), delarc_(), delnb_(), delnod_(), drwarc_(), edge_(), EMAN::EMData::extract_plane(), EMAN::EMData::extractline(), EMAN::Util::fftc_d(), fftc_d(), EMAN::Util::fftc_q(), fftc_q(), EMAN::Util::fftr_d(), fftr_d(), EMAN::Util::fftr_q(), fftr_q(), EMAN::EMData::find_3d_threshold(), EMAN::nnSSNR_ctfReconstructor::finish(), EMAN::nn4_ctfReconstructor::finish(), EMAN::nnSSNR_Reconstructor::finish(), EMAN::nn4Reconstructor::finish(), EMAN::fourierproduct(), EMAN::EMData::get_complex_at(), EMAN::TestUtil::get_pixel_value_by_dist1(), EMAN::TestUtil::get_pixel_value_by_dist2(), EMAN::Util::getBaldwinGridWeights(), EMAN::EMData::getconvpt2d_kbi0(), getnp_(), EMAN::Util::hypot_fast(), EMAN::Util::hypot_fast_int(), EMAN::FourierInserter3DMode5::insert_pixel(), EMAN::WienerFourierReconstructor::insert_slice(), EMAN::GaussFFTProjector::interp_ft_3d(), main(), EMAN::EMData::make_footprint(), max2d(), max3d(), EMAN::Util::multiref_polar_ali_2d_local(), EMAN::Util::multiref_polar_ali_2d_local_psi(), nearnd_(), EMAN::InterpolationFunctoidMode5::operate(), optim_(), EMAN::AutoMask3D2Processor::process_inplace(), EMAN::LinearPyramidProcessor::process_inplace(), EMAN::EMData::rotavg(), EMAN::EMData::rotavg_i(), EMAN::BoxingTools::set_region(), slamc2_(), swap_(), trfind_(), trlist_(), trmesh_(), EMAN::Util::trmsh3_(), trplot_(), and vrplot_().
| #define dabs | ( | x | ) | (doublereal)abs(x) |
Definition at line 47 of file lapackblas.h.
Referenced by isamax_(), sbdsqr_(), slae2_(), slaed0_(), slaed2_(), slaed4_(), slaed5_(), slaed6_(), slaed8_(), slaev2_(), slamc2_(), slange_(), slanst_(), slansy_(), slapy2_(), slarfg_(), slartg_(), slas2_(), slascl_(), slasq1_(), slasq2_(), slasq3_(), slassq_(), slasv2_(), snrm2_(), sstedc_(), ssteqr_(), and ssterf_().
| #define df2cmax | ( | a, | |||
| b | ) | (doublereal)f2cmax(a,b) |
| #define df2cmin | ( | a, | |||
| b | ) | (doublereal)f2cmin(a,b) |
| #define f2cmax | ( | a, | |||
| b | ) | ((a) >= (b) ? (a) : (b)) |
Definition at line 49 of file lapackblas.h.
Referenced by sbdsqr_(), sgebd2_(), sgebrd_(), sgelq2_(), sgelqf_(), sgemm_(), sgemv_(), sgeqr2_(), sgeqrf_(), sger_(), sgesvd_(), slaed0_(), slaed1_(), slaed2_(), slaed3_(), slaed6_(), slaed7_(), slaed8_(), slaed9_(), slamc2_(), slascl_(), slasq2_(), slasq3_(), slasr_(), sorg2l_(), sorg2r_(), sorgbr_(), sorgl2_(), sorglq_(), sorgql_(), sorgqr_(), sorgtr_(), sorm2r_(), sormbr_(), sorml2_(), sormlq_(), sormqr_(), sstedc_(), ssteqr_(), ssyev_(), ssymv_(), ssyr2_(), ssyr2k_(), ssytd2_(), ssytrd_(), strmm_(), and strmv_().
| #define f2cmin | ( | a, | |||
| b | ) | ((a) <= (b) ? (a) : (b)) |
Definition at line 48 of file lapackblas.h.
Referenced by ilaenv_(), sgebd2_(), sgebrd_(), sgelq2_(), sgelqf_(), sgeqr2_(), sgeqrf_(), sgesvd_(), slabrd_(), slacpy_(), slaed1_(), slaed2_(), slaed7_(), slaed8_(), slamc2_(), slange_(), slascl_(), slaset_(), slasq3_(), slatrd_(), sorgbr_(), sorglq_(), sorgql_(), sorgqr_(), sormbr_(), sormlq_(), sormqr_(), and ssytd2_().
| #define FALSE_ (0) |
| #define TRUE_ (1) |
Typedef Documentation
| typedef double doublereal |
Definition at line 36 of file lapackblas.h.
| typedef short flag |
Definition at line 39 of file lapackblas.h.
| typedef short ftnint |
Definition at line 41 of file lapackblas.h.
| typedef short ftnlen |
Definition at line 40 of file lapackblas.h.
| typedef int integer |
Definition at line 34 of file lapackblas.h.
| typedef int logical |
Definition at line 37 of file lapackblas.h.
| typedef float real |
Definition at line 35 of file lapackblas.h.
Function Documentation
| real drand | ( | void | ) |
Definition at line 56 of file lapackblas.cpp.
References integer.
Referenced by ilaenv_().
00057 { 00058 /* -- LAPACK auxiliary routine (version 3.0) -- 00059 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00060 Courant Institute, Argonne National Lab, and Rice University 00061 June 30, 1998 00062 00063 00064 Purpose 00065 ======= 00066 00067 IEEECK is called from the ILAENV to verify that Infinity and 00068 possibly NaN arithmetic is safe (i.e. will not trap). 00069 00070 Arguments 00071 ========= 00072 00073 ISPEC (input) INTEGER 00074 Specifies whether to test just for inifinity arithmetic 00075 or whether to test for infinity and NaN arithmetic. 00076 = 0: Verify infinity arithmetic only. 00077 = 1: Verify infinity and NaN arithmetic. 00078 00079 ZERO (input) REAL 00080 Must contain the value 0.0 00081 This is passed to prevent the compiler from optimizing 00082 away this code. 00083 00084 ONE (input) REAL 00085 Must contain the value 1.0 00086 This is passed to prevent the compiler from optimizing 00087 away this code. 00088 00089 RETURN VALUE: INTEGER 00090 = 0: Arithmetic failed to produce the correct answers 00091 = 1: Arithmetic produced the correct answers */ 00092 /* System generated locals */ 00093 integer ret_val; 00094 /* Local variables */ 00095 static real neginf, posinf, negzro, newzro, nan1, nan2, nan3, nan4, nan5, 00096 nan6; 00097 00098 00099 ret_val = 1; 00100 00101 posinf = *one / *zero; 00102 if (posinf <= *one) { 00103 ret_val = 0; 00104 return ret_val; 00105 } 00106 00107 neginf = -(*one) / *zero; 00108 if (neginf >= *zero) { 00109 ret_val = 0; 00110 return ret_val; 00111 } 00112 00113 negzro = *one / (neginf + *one); 00114 if (negzro != *zero) { 00115 ret_val = 0; 00116 return ret_val; 00117 } 00118 00119 neginf = *one / negzro; 00120 if (neginf >= *zero) { 00121 ret_val = 0; 00122 return ret_val; 00123 } 00124 00125 newzro = negzro + *zero; 00126 if (newzro != *zero) { 00127 ret_val = 0; 00128 return ret_val; 00129 } 00130 00131 posinf = *one / newzro; 00132 if (posinf <= *one) { 00133 ret_val = 0; 00134 return ret_val; 00135 } 00136 00137 neginf *= posinf; 00138 if (neginf >= *zero) { 00139 ret_val = 0; 00140 return ret_val; 00141 } 00142 00143 posinf *= posinf; 00144 if (posinf <= *one) { 00145 ret_val = 0; 00146 return ret_val; 00147 } 00148 00149 00150 00151 00152 /* Return if we were only asked to check infinity arithmetic */ 00153 00154 if (*ispec == 0) { 00155 return ret_val; 00156 } 00157 00158 nan1 = posinf + neginf; 00159 00160 nan2 = posinf / neginf; 00161 00162 nan3 = posinf / posinf; 00163 00164 nan4 = posinf * *zero; 00165 00166 nan5 = neginf * negzro; 00167 00168 nan6 = nan5 * 0.f; 00169 00170 if (nan1 == nan1) { 00171 ret_val = 0; 00172 return ret_val; 00173 } 00174 00175 if (nan2 == nan2) { 00176 ret_val = 0; 00177 return ret_val; 00178 } 00179 00180 if (nan3 == nan3) { 00181 ret_val = 0; 00182 return ret_val; 00183 } 00184 00185 if (nan4 == nan4) { 00186 ret_val = 0; 00187 return ret_val; 00188 } 00189 00190 if (nan5 == nan5) { 00191 ret_val = 0; 00192 return ret_val; 00193 } 00194 00195 if (nan6 == nan6) { 00196 ret_val = 0; 00197 return ret_val; 00198 } 00199 00200 return ret_val; 00201 } /* ieeeck_ */
| integer ilaenv_ | ( | integer * | ispec, | |
| const char * | name__, | |||
| const char * | opts, | |||
| integer * | n1, | |||
| integer * | n2, | |||
| integer * | n3, | |||
| integer * | n4, | |||
| ftnlen | name_len, | |||
| ftnlen | opts_len | |||
| ) |
Definition at line 206 of file lapackblas.cpp.
References c__0, c__1, f2cmin, ieeeck_(), integer, nx, s_cmp(), and s_copy().
Referenced by sgebrd_(), sgelqf_(), sgeqrf_(), sgesvd_(), slaed0_(), slasq2_(), sorgbr_(), sorglq_(), sorgql_(), sorgqr_(), sorgtr_(), sormbr_(), sormlq_(), sormqr_(), sstedc_(), ssyev_(), and ssytrd_().
00208 { 00209 /* -- LAPACK auxiliary routine (version 3.0) -- 00210 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00211 Courant Institute, Argonne National Lab, and Rice University 00212 June 30, 1999 00213 00214 00215 Purpose 00216 ======= 00217 00218 ILAENV is called from the LAPACK routines to choose problem-dependent 00219 parameters for the local environment. See ISPEC for a description of 00220 the parameters. 00221 00222 This version provides a set of parameters which should give good, 00223 but not optimal, performance on many of the currently available 00224 computers. Users are encouraged to modify this subroutine to set 00225 the tuning parameters for their particular machine using the option 00226 and problem size information in the arguments. 00227 00228 This routine will not function correctly if it is converted to all 00229 lower case. Converting it to all upper case is allowed. 00230 00231 Arguments 00232 ========= 00233 00234 ISPEC (input) INTEGER 00235 Specifies the parameter to be returned as the value of 00236 ILAENV. 00237 = 1: the optimal blocksize; if this value is 1, an unblocked 00238 algorithm will give the best performance. 00239 = 2: the minimum block size for which the block routine 00240 should be used; if the usable block size is less than 00241 this value, an unblocked routine should be used. 00242 = 3: the crossover point (in a block routine, for N less 00243 than this value, an unblocked routine should be used) 00244 = 4: the number of shifts, used in the nonsymmetric 00245 eigenvalue routines 00246 = 5: the MINIMUM Column dimension for blocking to be used; 00247 rectangular blocks must have dimension at least k by m, 00248 where k is given by ILAENV(2,...) and m by ILAENV(5,...) 00249 = 6: the crossover point for the SVD (when reducing an m by n 00250 matrix to bidiagonal form, if f2cmax(m,n)/min(m,n) exceeds 00251 this value, a QR factorization is used first to reduce 00252 the matrix to a triangular form.) 00253 = 7: the number of processors 00254 = 8: the crossover point for the multishift QR and QZ methods 00255 for nonsymmetric eigenvalue problems. 00256 = 9: maximum size of the subproblems at the bottom of the 00257 computation tree in the divide-and-conquer algorithm 00258 (used by xGELSD and xGESDD) 00259 =10: ieee NaN arithmetic can be trusted not to trap 00260 =11: infinity arithmetic can be trusted not to trap 00261 00262 NAME (input) CHARACTER*(*) 00263 The name of the calling subroutine, in either upper case or 00264 lower case. 00265 00266 OPTS (input) CHARACTER*(*) 00267 The character options to the subroutine NAME, concatenated 00268 into a single character string. For example, UPLO = 'U', 00269 TRANS = 'T', and DIAG = 'N' for a triangular routine would 00270 be specified as OPTS = 'UTN'. 00271 00272 N1 (input) INTEGER 00273 N2 (input) INTEGER 00274 N3 (input) INTEGER 00275 N4 (input) INTEGER 00276 Problem dimensions for the subroutine NAME; these may not all 00277 be required. 00278 00279 (ILAENV) (output) INTEGER 00280 >= 0: the value of the parameter specified by ISPEC 00281 < 0: if ILAENV = -k, the k-th argument had an illegal value. 00282 00283 Further Details 00284 =============== 00285 00286 The following conventions have been used when calling ILAENV from the 00287 LAPACK routines: 00288 1) OPTS is a concatenation of all of the character options to 00289 subroutine NAME, in the same order that they appear in the 00290 argument list for NAME, even if they are not used in determining 00291 the value of the parameter specified by ISPEC. 00292 2) The problem dimensions N1, N2, N3, N4 are specified in the order 00293 that they appear in the argument list for NAME. N1 is used 00294 first, N2 second, and so on, and unused problem dimensions are 00295 passed a value of -1. 00296 3) The parameter value returned by ILAENV is checked for validity in 00297 the calling subroutine. For example, ILAENV is used to retrieve 00298 the optimal blocksize for STRTRI as follows: 00299 00300 NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) 00301 IF( NB.LE.1 ) NB = MAX( 1, N ) 00302 00303 ===================================================================== */ 00304 /* Table of constant values */ 00305 static integer c__0 = 0; 00306 static real c_b162 = 0.f; 00307 static real c_b163 = 1.f; 00308 static integer c__1 = 1; 00309 00310 /* System generated locals */ 00311 integer ret_val; 00312 /* Builtin functions 00313 Subroutine */ void s_copy(char *, const char *, ftnlen, ftnlen); 00314 integer s_cmp(char *, const char *, ftnlen, ftnlen); 00315 /* Local variables */ 00316 static integer i__; 00317 static logical cname, sname; 00318 static integer nbmin; 00319 static char c1[1], c2[2], c3[3], c4[2]; 00320 static integer ic, nb; 00321 extern integer ieeeck_(integer *, real *, real *); 00322 static integer iz, nx; 00323 static char subnam[6]; 00324 00325 00326 00327 00328 switch (*ispec) { 00329 case 1: goto L100; 00330 case 2: goto L100; 00331 case 3: goto L100; 00332 case 4: goto L400; 00333 case 5: goto L500; 00334 case 6: goto L600; 00335 case 7: goto L700; 00336 case 8: goto L800; 00337 case 9: goto L900; 00338 case 10: goto L1000; 00339 case 11: goto L1100; 00340 } 00341 00342 /* Invalid value for ISPEC */ 00343 00344 ret_val = -1; 00345 return ret_val; 00346 00347 L100: 00348 00349 /* Convert NAME to upper case if the first character is lower case. */ 00350 00351 ret_val = 1; 00352 s_copy(subnam, name__, (ftnlen)6, name_len); 00353 ic = *(unsigned char *)subnam; 00354 iz = 'Z'; 00355 if (iz == 90 || iz == 122) { 00356 00357 /* ASCII character set */ 00358 00359 if (ic >= 97 && ic <= 122) { 00360 *(unsigned char *)subnam = (char) (ic - 32); 00361 for (i__ = 2; i__ <= 6; ++i__) { 00362 ic = *(unsigned char *)&subnam[i__ - 1]; 00363 if (ic >= 97 && ic <= 122) { 00364 *(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32); 00365 } 00366 /* L10: */ 00367 } 00368 } 00369 00370 } else if (iz == 233 || iz == 169) { 00371 00372 /* EBCDIC character set */ 00373 00374 if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 || ic >= 162 && 00375 ic <= 169) { 00376 *(unsigned char *)subnam = (char) (ic + 64); 00377 for (i__ = 2; i__ <= 6; ++i__) { 00378 ic = *(unsigned char *)&subnam[i__ - 1]; 00379 if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 || ic >= 00380 162 && ic <= 169) { 00381 *(unsigned char *)&subnam[i__ - 1] = (char) (ic + 64); 00382 } 00383 /* L20: */ 00384 } 00385 } 00386 00387 } else if (iz == 218 || iz == 250) { 00388 00389 /* Prime machines: ASCII+128 */ 00390 00391 if (ic >= 225 && ic <= 250) { 00392 *(unsigned char *)subnam = (char) (ic - 32); 00393 for (i__ = 2; i__ <= 6; ++i__) { 00394 ic = *(unsigned char *)&subnam[i__ - 1]; 00395 if (ic >= 225 && ic <= 250) { 00396 *(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32); 00397 } 00398 /* L30: */ 00399 } 00400 } 00401 } 00402 00403 *(unsigned char *)c1 = *(unsigned char *)subnam; 00404 sname = *(unsigned char *)c1 == 'S' || *(unsigned char *)c1 == 'D'; 00405 cname = *(unsigned char *)c1 == 'C' || *(unsigned char *)c1 == 'Z'; 00406 if (! (cname || sname)) { 00407 return ret_val; 00408 } 00409 s_copy(c2, subnam + 1, (ftnlen)2, (ftnlen)2); 00410 s_copy(c3, subnam + 3, (ftnlen)3, (ftnlen)3); 00411 s_copy(c4, c3 + 1, (ftnlen)2, (ftnlen)2); 00412 00413 switch (*ispec) { 00414 case 1: goto L110; 00415 case 2: goto L200; 00416 case 3: goto L300; 00417 } 00418 00419 L110: 00420 00421 /* ISPEC = 1: block size 00422 00423 In these examples, separate code is provided for setting NB for 00424 real and complex. We assume that NB will take the same value in 00425 single or double precision. */ 00426 00427 nb = 1; 00428 00429 if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) { 00430 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00431 if (sname) { 00432 nb = 64; 00433 } else { 00434 nb = 64; 00435 } 00436 } else if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, 00437 "RQF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen) 00438 3, (ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3) 00439 == 0) { 00440 if (sname) { 00441 nb = 32; 00442 } else { 00443 nb = 32; 00444 } 00445 } else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) { 00446 if (sname) { 00447 nb = 32; 00448 } else { 00449 nb = 32; 00450 } 00451 } else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) { 00452 if (sname) { 00453 nb = 32; 00454 } else { 00455 nb = 32; 00456 } 00457 } else if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) { 00458 if (sname) { 00459 nb = 64; 00460 } else { 00461 nb = 64; 00462 } 00463 } 00464 } else if (s_cmp(c2, "PO", (ftnlen)2, (ftnlen)2) == 0) { 00465 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00466 if (sname) { 00467 nb = 64; 00468 } else { 00469 nb = 64; 00470 } 00471 } 00472 } else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) { 00473 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00474 if (sname) { 00475 nb = 64; 00476 } else { 00477 nb = 64; 00478 } 00479 } else if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) { 00480 nb = 32; 00481 } else if (sname && s_cmp(c3, "GST", (ftnlen)3, (ftnlen)3) == 0) { 00482 nb = 64; 00483 } 00484 } else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) { 00485 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00486 nb = 64; 00487 } else if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) { 00488 nb = 32; 00489 } else if (s_cmp(c3, "GST", (ftnlen)3, (ftnlen)3) == 0) { 00490 nb = 64; 00491 } 00492 } else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) { 00493 if (*(unsigned char *)c3 == 'G') { 00494 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00495 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00496 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00497 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00498 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00499 ftnlen)2, (ftnlen)2) == 0) { 00500 nb = 32; 00501 } 00502 } else if (*(unsigned char *)c3 == 'M') { 00503 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00504 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00505 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00506 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00507 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00508 ftnlen)2, (ftnlen)2) == 0) { 00509 nb = 32; 00510 } 00511 } 00512 } else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) { 00513 if (*(unsigned char *)c3 == 'G') { 00514 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00515 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00516 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00517 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00518 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00519 ftnlen)2, (ftnlen)2) == 0) { 00520 nb = 32; 00521 } 00522 } else if (*(unsigned char *)c3 == 'M') { 00523 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00524 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00525 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00526 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00527 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00528 ftnlen)2, (ftnlen)2) == 0) { 00529 nb = 32; 00530 } 00531 } 00532 } else if (s_cmp(c2, "GB", (ftnlen)2, (ftnlen)2) == 0) { 00533 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00534 if (sname) { 00535 if (*n4 <= 64) { 00536 nb = 1; 00537 } else { 00538 nb = 32; 00539 } 00540 } else { 00541 if (*n4 <= 64) { 00542 nb = 1; 00543 } else { 00544 nb = 32; 00545 } 00546 } 00547 } 00548 } else if (s_cmp(c2, "PB", (ftnlen)2, (ftnlen)2) == 0) { 00549 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00550 if (sname) { 00551 if (*n2 <= 64) { 00552 nb = 1; 00553 } else { 00554 nb = 32; 00555 } 00556 } else { 00557 if (*n2 <= 64) { 00558 nb = 1; 00559 } else { 00560 nb = 32; 00561 } 00562 } 00563 } 00564 } else if (s_cmp(c2, "TR", (ftnlen)2, (ftnlen)2) == 0) { 00565 if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) { 00566 if (sname) { 00567 nb = 64; 00568 } else { 00569 nb = 64; 00570 } 00571 } 00572 } else if (s_cmp(c2, "LA", (ftnlen)2, (ftnlen)2) == 0) { 00573 if (s_cmp(c3, "UUM", (ftnlen)3, (ftnlen)3) == 0) { 00574 if (sname) { 00575 nb = 64; 00576 } else { 00577 nb = 64; 00578 } 00579 } 00580 } else if (sname && s_cmp(c2, "ST", (ftnlen)2, (ftnlen)2) == 0) { 00581 if (s_cmp(c3, "EBZ", (ftnlen)3, (ftnlen)3) == 0) { 00582 nb = 1; 00583 } 00584 } 00585 ret_val = nb; 00586 return ret_val; 00587 00588 L200: 00589 00590 /* ISPEC = 2: minimum block size */ 00591 00592 nbmin = 2; 00593 if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) { 00594 if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "RQF", ( 00595 ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)3, ( 00596 ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3) == 0) 00597 { 00598 if (sname) { 00599 nbmin = 2; 00600 } else { 00601 nbmin = 2; 00602 } 00603 } else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) { 00604 if (sname) { 00605 nbmin = 2; 00606 } else { 00607 nbmin = 2; 00608 } 00609 } else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) { 00610 if (sname) { 00611 nbmin = 2; 00612 } else { 00613 nbmin = 2; 00614 } 00615 } else if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) { 00616 if (sname) { 00617 nbmin = 2; 00618 } else { 00619 nbmin = 2; 00620 } 00621 } 00622 } else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) { 00623 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) { 00624 if (sname) { 00625 nbmin = 8; 00626 } else { 00627 nbmin = 8; 00628 } 00629 } else if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) { 00630 nbmin = 2; 00631 } 00632 } else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) { 00633 if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) { 00634 nbmin = 2; 00635 } 00636 } else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) { 00637 if (*(unsigned char *)c3 == 'G') { 00638 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00639 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00640 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00641 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00642 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00643 ftnlen)2, (ftnlen)2) == 0) { 00644 nbmin = 2; 00645 } 00646 } else if (*(unsigned char *)c3 == 'M') { 00647 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00648 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00649 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00650 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00651 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00652 ftnlen)2, (ftnlen)2) == 0) { 00653 nbmin = 2; 00654 } 00655 } 00656 } else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) { 00657 if (*(unsigned char *)c3 == 'G') { 00658 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00659 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00660 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00661 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00662 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00663 ftnlen)2, (ftnlen)2) == 0) { 00664 nbmin = 2; 00665 } 00666 } else if (*(unsigned char *)c3 == 'M') { 00667 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00668 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00669 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00670 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00671 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00672 ftnlen)2, (ftnlen)2) == 0) { 00673 nbmin = 2; 00674 } 00675 } 00676 } 00677 ret_val = nbmin; 00678 return ret_val; 00679 00680 L300: 00681 00682 /* ISPEC = 3: crossover point */ 00683 00684 nx = 0; 00685 if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) { 00686 if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "RQF", ( 00687 ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)3, ( 00688 ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3) == 0) 00689 { 00690 if (sname) { 00691 nx = 128; 00692 } else { 00693 nx = 128; 00694 } 00695 } else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) { 00696 if (sname) { 00697 nx = 128; 00698 } else { 00699 nx = 128; 00700 } 00701 } else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) { 00702 if (sname) { 00703 nx = 128; 00704 } else { 00705 nx = 128; 00706 } 00707 } 00708 } else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) { 00709 if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) { 00710 nx = 32; 00711 } 00712 } else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) { 00713 if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) { 00714 nx = 32; 00715 } 00716 } else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) { 00717 if (*(unsigned char *)c3 == 'G') { 00718 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00719 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00720 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00721 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00722 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00723 ftnlen)2, (ftnlen)2) == 0) { 00724 nx = 128; 00725 } 00726 } 00727 } else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) { 00728 if (*(unsigned char *)c3 == 'G') { 00729 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ", 00730 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, ( 00731 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) == 00732 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp( 00733 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", ( 00734 ftnlen)2, (ftnlen)2) == 0) { 00735 nx = 128; 00736 } 00737 } 00738 } 00739 ret_val = nx; 00740 return ret_val; 00741 00742 L400: 00743 00744 /* ISPEC = 4: number of shifts (used by xHSEQR) */ 00745 00746 ret_val = 6; 00747 return ret_val; 00748 00749 L500: 00750 00751 /* ISPEC = 5: minimum column dimension (not used) */ 00752 00753 ret_val = 2; 00754 return ret_val; 00755 00756 L600: 00757 00758 /* ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) */ 00759 00760 ret_val = (integer) ((real) f2cmin(*n1,*n2) * 1.6f); 00761 return ret_val; 00762 00763 L700: 00764 00765 /* ISPEC = 7: number of processors (not used) */ 00766 00767 ret_val = 1; 00768 return ret_val; 00769 00770 L800: 00771 00772 /* ISPEC = 8: crossover point for multishift (used by xHSEQR) */ 00773 00774 ret_val = 50; 00775 return ret_val; 00776 00777 L900: 00778 00779 /* ISPEC = 9: maximum size of the subproblems at the bottom of the 00780 computation tree in the divide-and-conquer algorithm 00781 (used by xGELSD and xGESDD) */ 00782 00783 ret_val = 25; 00784 return ret_val; 00785 00786 L1000: 00787 00788 /* ISPEC = 10: ieee NaN arithmetic can be trusted not to trap 00789 00790 ILAENV = 0 */ 00791 ret_val = 1; 00792 if (ret_val == 1) { 00793 ret_val = ieeeck_(&c__0, &c_b162, &c_b163); 00794 } 00795 return ret_val; 00796 00797 L1100: 00798 00799 /* ISPEC = 11: infinity arithmetic can be trusted not to trap 00800 00801 ILAENV = 0 */ 00802 ret_val = 1; 00803 if (ret_val == 1) { 00804 ret_val = ieeeck_(&c__1, &c_b162, &c_b163); 00805 } 00806 return ret_val; 00807 00808 /* End of ILAENV */ 00809 00810 } /* ilaenv_ */
| logical lsame_ | ( | const char * | ca, | |
| const char * | cb | |||
| ) |
Definition at line 814 of file lapackblas.cpp.
References integer.
Referenced by sbdsqr_(), sgemm_(), sgemv_(), sgesvd_(), slacpy_(), slamch_(), slange_(), slanst_(), slansy_(), slarf_(), slarfb_(), slarft_(), slascl_(), slaset_(), slasr_(), slasrt_(), slatrd_(), sorgbr_(), sorgtr_(), sorm2r_(), sormbr_(), sorml2_(), sormlq_(), sormqr_(), sstedc_(), ssteqr_(), sstevd_(), ssyev_(), ssymv_(), ssyr2_(), ssyr2k_(), ssytd2_(), ssytrd_(), strmm_(), and strmv_().
00815 { 00816 /* -- LAPACK auxiliary routine (version 3.0) -- 00817 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00818 Courant Institute, Argonne National Lab, and Rice University 00819 September 30, 1994 00820 00821 00822 Purpose 00823 ======= 00824 00825 LSAME returns .TRUE. if CA is the same letter as CB regardless of 00826 case. 00827 00828 Arguments 00829 ========= 00830 00831 CA (input) CHARACTER*1 00832 CB (input) CHARACTER*1 00833 CA and CB specify the single characters to be compared. 00834 00835 ===================================================================== 00836 00837 00838 00839 Test if the characters are equal */ 00840 /* System generated locals */ 00841 logical ret_val; 00842 /* Local variables */ 00843 static integer inta, intb, zcode; 00844 00845 00846 ret_val = *(unsigned char *)ca == *(unsigned char *)cb; 00847 if (ret_val) { 00848 return ret_val; 00849 } 00850 00851 /* Now test for equivalence if both characters are alphabetic. */ 00852 00853 zcode = 'Z'; 00854 00855 /* Use 'Z' rather than 'A' so that ASCII can be detected on Prime 00856 machines, on which ICHAR returns a value with bit 8 set. 00857 ICHAR('A') on Prime machines returns 193 which is the same as 00858 ICHAR('A') on an EBCDIC machine. */ 00859 00860 inta = *(unsigned char *)ca; 00861 intb = *(unsigned char *)cb; 00862 00863 if (zcode == 90 || zcode == 122) { 00864 00865 /* ASCII is assumed - ZCODE is the ASCII code of either lower o 00866 r 00867 upper case 'Z'. */ 00868 00869 if (inta >= 97 && inta <= 122) { 00870 inta += -32; 00871 } 00872 if (intb >= 97 && intb <= 122) { 00873 intb += -32; 00874 } 00875 00876 } else if (zcode == 233 || zcode == 169) { 00877 00878 /* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower 00879 or 00880 upper case 'Z'. */ 00881 00882 if (inta >= 129 && inta <= 137 || inta >= 145 && inta <= 153 || inta 00883 >= 162 && inta <= 169) { 00884 inta += 64; 00885 } 00886 if (intb >= 129 && intb <= 137 || intb >= 145 && intb <= 153 || intb 00887 >= 162 && intb <= 169) { 00888 intb += 64; 00889 } 00890 00891 } else if (zcode == 218 || zcode == 250) { 00892 00893 /* ASCII is assumed, on Prime machines - ZCODE is the ASCII cod 00894 e 00895 plus 128 of either lower or upper case 'Z'. */ 00896 00897 if (inta >= 225 && inta <= 250) { 00898 inta += -32; 00899 } 00900 if (intb >= 225 && intb <= 250) { 00901 intb += -32; 00902 } 00903 } 00904 ret_val = inta == intb; 00905 00906 /* RETURN 00907 00908 End of LSAME */ 00909 00910 return ret_val; 00911 } /* lsame_ */
| double pow_ri | ( | real * | ap, | |
| integer * | bp | |||
| ) |
Definition at line 918 of file lapackblas.cpp.
Referenced by slaed6_(), slamc2_(), slamch_(), and slartg_().
00920 { 00921 double pow, x; 00922 integer n; 00923 unsigned long u; 00924 00925 pow = 1; 00926 x = *ap; 00927 n = *bp; 00928 00929 if(n != 0) 00930 { 00931 if(n < 0) 00932 { 00933 n = -n; 00934 x = 1/x; 00935 } 00936 for(u = n; ; ) 00937 { 00938 if(u & 01) 00939 pow *= x; 00940 if(u >>= 1) 00941 x *= x; 00942 else 00943 break; 00944 } 00945 } 00946 return(pow); 00947 }
| double r_sign | ( | real * | a, | |
| real * | b | |||
| ) |
Definition at line 37 of file lapackblas.cpp.
Referenced by sgesvd_(), sormbr_(), sormlq_(), and sormqr_().
00039 { 00040 ftnlen i, n, nc; 00041 const char *f__rp; 00042 00043 n = (int)*np; 00044 for(i = 0 ; i < n ; ++i) { 00045 nc = ll; 00046 if(rnp[i] < nc) nc = rnp[i]; 00047 ll -= nc; 00048 f__rp = rpp[i]; 00049 while(--nc >= 0) *lp++ = *f__rp++; 00050 } 00051 while(--ll >= 0) 00052 *lp++ = ' '; 00053 return 0; 00054 }
Definition at line 1071 of file lapackblas.cpp.
Referenced by dcsrch_(), ilaenv_(), lnsrlb_(), mainlb_(), prn3lb_(), and setulb_().
01073 { 01074 register unsigned char *a, *aend, *b, *bend; 01075 a = (unsigned char *)a0; 01076 b = (unsigned char *)b0; 01077 aend = a + la; 01078 bend = b + lb; 01079 01080 if(la <= lb) 01081 { 01082 while(a < aend) 01083 if(*a != *b) 01084 return( *a - *b ); 01085 else 01086 { ++a; ++b; } 01087 01088 while(b < bend) 01089 if(*b != ' ') 01090 return( ' ' - *b ); 01091 else ++b; 01092 } 01093 01094 else 01095 { 01096 while(b < bend) 01097 if(*a == *b) 01098 { ++a; ++b; } 01099 else 01100 return( *a - *b ); 01101 while(a < aend) 01102 if(*a != ' ') 01103 return(*a - ' '); 01104 else ++a; 01105 } 01106 return(0); 01107 }
Definition at line 1121 of file lapackblas.cpp.
Referenced by dcsrch_(), errclb_(), ilaenv_(), lnsrlb_(), mainlb_(), and prn2lb_().
01123 { 01124 register char *aend; 01125 const register char *bend; 01126 01127 aend = a + la; 01128 01129 if(la <= lb) 01130 #ifndef NO_OVERWRITE 01131 if (a <= b || a >= b + la) 01132 #endif 01133 while(a < aend) 01134 *a++ = *b++; 01135 #ifndef NO_OVERWRITE 01136 else 01137 for(b += la; a < aend; ) 01138 *--aend = *--b; 01139 #endif 01140 01141 else { 01142 bend = b + lb; 01143 #ifndef NO_OVERWRITE 01144 if (a <= b || a >= bend) 01145 #endif 01146 while(b < bend) 01147 *a++ = *b++; 01148 #ifndef NO_OVERWRITE 01149 else { 01150 a += lb; 01151 while(b < bend) 01152 *--a = *--bend; 01153 a += lb; 01154 } 01155 #endif 01156 while(a < aend) 01157 *a++ = ' '; 01158 } 01159 }
Definition at line 994 of file lapackblas.cpp.
References integer.
Referenced by EMAN::PCA::dopca_ooc(), EMAN::PCA::Lanczos(), EMAN::PCAlarge::Lanczos(), EMAN::PCA::Lanczos_ooc(), slatrd_(), and ssytd2_().
00996 { 00997 /* System generated locals */ 00998 integer i__1; 00999 /* Local variables */ 01000 static integer i__, m, ix, iy, mp1; 01001 /* constant times a vector plus a vector. 01002 uses unrolled loop for increments equal to one. 01003 jack dongarra, linpack, 3/11/78. 01004 modified 12/3/93, array(1) declarations changed to array(*) 01005 Parameter adjustments */ 01006 --sy; 01007 --sx; 01008 /* Function Body */ 01009 if (*n <= 0) { 01010 return 0; 01011 } 01012 if (*sa == 0.f) { 01013 return 0; 01014 } 01015 if (*incx == 1 && *incy == 1) { 01016 goto L20; 01017 } 01018 /* code for unequal increments or equal increments 01019 not equal to 1 */ 01020 ix = 1; 01021 iy = 1; 01022 if (*incx < 0) { 01023 ix = (-(*n) + 1) * *incx + 1; 01024 } 01025 if (*incy < 0) { 01026 iy = (-(*n) + 1) * *incy + 1; 01027 } 01028 i__1 = *n; 01029 for (i__ = 1; i__ <= i__1; ++i__) { 01030 sy[iy] += *sa * sx[ix]; 01031 ix += *incx; 01032 iy += *incy; 01033 /* L10: */ 01034 } 01035 return 0; 01036 /* code for both increments equal to 1 01037 clean-up loop */ 01038 L20: 01039 m = *n % 4; 01040 if (m == 0) { 01041 goto L40; 01042 } 01043 i__1 = m; 01044 for (i__ = 1; i__ <= i__1; ++i__) { 01045 sy[i__] += *sa * sx[i__]; 01046 /* L30: */ 01047 } 01048 if (*n < 4) { 01049 return 0; 01050 } 01051 L40: 01052 mp1 = m + 1; 01053 i__1 = *n; 01054 for (i__ = mp1; i__ <= i__1; i__ += 4) { 01055 sy[i__] += *sa * sx[i__]; 01056 sy[i__ + 1] += *sa * sx[i__ + 1]; 01057 sy[i__ + 2] += *sa * sx[i__ + 2]; 01058 sy[i__ + 3] += *sa * sx[i__ + 3]; 01059 /* L50: */ 01060 } 01061 return 0; 01062 } /* saxpy_ */
| int sbdsqr_ | ( | const char * | uplo, | |
| integer * | n, | |||
| integer * | ncvt, | |||
| integer * | nru, | |||
| integer * | ncc, | |||
| real * | d__, | |||
| real * | e, | |||
| real * | vt, | |||
| integer * | ldvt, | |||
| real * | u, | |||
| integer * | ldu, | |||
| real * | c__, | |||
| integer * | ldc, | |||
| real * | work, | |||
| integer * | info | |||
| ) |
Definition at line 22954 of file lapackblas.cpp.
References c___ref, dabs, df2cmax, df2cmin, f2cmax, integer, lsame_(), r_sign(), slamch_(), slartg_(), slas2_(), slasq1_(), slasr_(), slasv2_(), sqrt(), srot_(), sscal_(), sswap_(), u_ref, vt_ref, and xerbla_().
Referenced by sgesvd_().
22957 { 22958 /* System generated locals */ 22959 integer c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1, 22960 i__2; 22961 real r__1, r__2, r__3, r__4; 22962 doublereal d__1; 22963 22964 /* Builtin functions */ 22965 // double pow_dd(doublereal *, doublereal *), sqrt(doublereal), r_sign(real *, real *); 22966 22967 /* Local variables */ 22968 static real abse; 22969 static integer idir; 22970 static real abss; 22971 static integer oldm; 22972 static real cosl; 22973 static integer isub, iter; 22974 static real unfl, sinl, cosr, smin, smax, sinr; 22975 extern /* Subroutine */ int srot_(integer *, real *, integer *, real *, 22976 integer *, real *, real *), slas2_(real *, real *, real *, real *, 22977 real *); 22978 static real f, g, h__; 22979 static integer i__, j, m; 22980 static real r__; 22981 extern logical lsame_(const char *, const char *); 22982 static real oldcs; 22983 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *); 22984 static integer oldll; 22985 static real shift, sigmn, oldsn; 22986 static integer maxit; 22987 static real sminl; 22988 extern /* Subroutine */ int slasr_(const char *, const char *, const char *, integer *, 22989 integer *, real *, real *, real *, integer *); 22990 static real sigmx; 22991 static logical lower; 22992 extern /* Subroutine */ int sswap_(integer *, real *, integer *, real *, 22993 integer *), slasq1_(integer *, real *, real *, real *, integer *), 22994 slasv2_(real *, real *, real *, real *, real *, real *, real *, 22995 real *, real *); 22996 static real cs; 22997 static integer ll; 22998 static real sn, mu; 22999 extern doublereal slamch_(const char *); 23000 extern /* Subroutine */ int xerbla_(const char *, integer *); 23001 static real sminoa; 23002 extern /* Subroutine */ int slartg_(real *, real *, real *, real *, real * 23003 ); 23004 static real thresh; 23005 static logical rotate; 23006 static real sminlo; 23007 static integer nm1; 23008 static real tolmul; 23009 static integer nm12, nm13, lll; 23010 static real eps, sll, tol; 23011 23012 23013 #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1] 23014 #define u_ref(a_1,a_2) u[(a_2)*u_dim1 + a_1] 23015 #define vt_ref(a_1,a_2) vt[(a_2)*vt_dim1 + a_1] 23016 23017 23018 /* -- LAPACK routine (version 3.0) -- 23019 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 23020 Courant Institute, Argonne National Lab, and Rice University 23021 October 31, 1999 23022 23023 23024 Purpose 23025 ======= 23026 23027 SBDSQR computes the singular value decomposition (SVD) of a real 23028 N-by-N (upper or lower) bidiagonal matrix B: B = Q * S * P' (P' 23029 denotes the transpose of P), where S is a diagonal matrix with 23030 non-negative diagonal elements (the singular values of B), and Q 23031 and P are orthogonal matrices. 23032 23033 The routine computes S, and optionally computes U * Q, P' * VT, 23034 or Q' * C, for given real input matrices U, VT, and C. 23035 23036 See "Computing Small Singular Values of Bidiagonal Matrices With 23037 Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, 23038 LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, 23039 no. 5, pp. 873-912, Sept 1990) and 23040 "Accurate singular values and differential qd algorithms," by 23041 B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics 23042 Department, University of California at Berkeley, July 1992 23043 for a detailed description of the algorithm. 23044 23045 Arguments 23046 ========= 23047 23048 UPLO (input) CHARACTER*1 23049 = 'U': B is upper bidiagonal; 23050 = 'L': B is lower bidiagonal. 23051 23052 N (input) INTEGER 23053 The order of the matrix B. N >= 0. 23054 23055 NCVT (input) INTEGER 23056 The number of columns of the matrix VT. NCVT >= 0. 23057 23058 NRU (input) INTEGER 23059 The number of rows of the matrix U. NRU >= 0. 23060 23061 NCC (input) INTEGER 23062 The number of columns of the matrix C. NCC >= 0. 23063 23064 D (input/output) REAL array, dimension (N) 23065 On entry, the n diagonal elements of the bidiagonal matrix B. 23066 On exit, if INFO=0, the singular values of B in decreasing 23067 order. 23068 23069 E (input/output) REAL array, dimension (N) 23070 On entry, the elements of E contain the 23071 offdiagonal elements of the bidiagonal matrix whose SVD 23072 is desired. On normal exit (INFO = 0), E is destroyed. 23073 If the algorithm does not converge (INFO > 0), D and E 23074 will contain the diagonal and superdiagonal elements of a 23075 bidiagonal matrix orthogonally equivalent to the one given 23076 as input. E(N) is used for workspace. 23077 23078 VT (input/output) REAL array, dimension (LDVT, NCVT) 23079 On entry, an N-by-NCVT matrix VT. 23080 On exit, VT is overwritten by P' * VT. 23081 VT is not referenced if NCVT = 0. 23082 23083 LDVT (input) INTEGER 23084 The leading dimension of the array VT. 23085 LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0. 23086 23087 U (input/output) REAL array, dimension (LDU, N) 23088 On entry, an NRU-by-N matrix U. 23089 On exit, U is overwritten by U * Q. 23090 U is not referenced if NRU = 0. 23091 23092 LDU (input) INTEGER 23093 The leading dimension of the array U. LDU >= max(1,NRU). 23094 23095 C (input/output) REAL array, dimension (LDC, NCC) 23096 On entry, an N-by-NCC matrix C. 23097 On exit, C is overwritten by Q' * C. 23098 C is not referenced if NCC = 0. 23099 23100 LDC (input) INTEGER 23101 The leading dimension of the array C. 23102 LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0. 23103 23104 WORK (workspace) REAL array, dimension (4*N) 23105 23106 INFO (output) INTEGER 23107 = 0: successful exit 23108 < 0: If INFO = -i, the i-th argument had an illegal value 23109 > 0: the algorithm did not converge; D and E contain the 23110 elements of a bidiagonal matrix which is orthogonally 23111 similar to the input matrix B; if INFO = i, i 23112 elements of E have not converged to zero. 23113 23114 Internal Parameters 23115 =================== 23116 23117 TOLMUL REAL, default = max(10,min(100,EPS**(-1/8))) 23118 TOLMUL controls the convergence criterion of the QR loop. 23119 If it is positive, TOLMUL*EPS is the desired relative 23120 precision in the computed singular values. 23121 If it is negative, abs(TOLMUL*EPS*sigma_max) is the 23122 desired absolute accuracy in the computed singular 23123 values (corresponds to relative accuracy 23124 abs(TOLMUL*EPS) in the largest singular value. 23125 abs(TOLMUL) should be between 1 and 1/EPS, and preferably 23126 between 10 (for fast convergence) and .1/EPS 23127 (for there to be some accuracy in the results). 23128 Default is to lose at either one eighth or 2 of the 23129 available decimal digits in each computed singular value 23130 (whichever is smaller). 23131 23132 MAXITR INTEGER, default = 6 23133 MAXITR controls the maximum number of passes of the 23134 algorithm through its inner loop. The algorithms stops 23135 (and so fails to converge) if the number of passes 23136 through the inner loop exceeds MAXITR*N**2. 23137 23138 ===================================================================== 23139 23140 23141 Test the input parameters. 23142 23143 Parameter adjustments */ 23144 --d__; 23145 --e; 23146 vt_dim1 = *ldvt; 23147 vt_offset = 1 + vt_dim1 * 1; 23148 vt -= vt_offset; 23149 u_dim1 = *ldu; 23150 u_offset = 1 + u_dim1 * 1; 23151 u -= u_offset; 23152 c_dim1 = *ldc; 23153 c_offset = 1 + c_dim1 * 1; 23154 c__ -= c_offset; 23155 --work; 23156 23157 /* Function Body */ 23158 *info = 0; 23159 lower = lsame_(uplo, "L"); 23160 if (! lsame_(uplo, "U") && ! lower) { 23161 *info = -1; 23162 } else if (*n < 0) { 23163 *info = -2; 23164 } else if (*ncvt < 0) { 23165 *info = -3; 23166 } else if (*nru < 0) { 23167 *info = -4; 23168 } else if (*ncc < 0) { 23169 *info = -5; 23170 } else if (*ncvt == 0 && *ldvt < 1 || *ncvt > 0 && *ldvt < f2cmax(1,*n)) { 23171 *info = -9; 23172 } else if (*ldu < f2cmax(1,*nru)) { 23173 *info = -11; 23174 } else if (*ncc == 0 && *ldc < 1 || *ncc > 0 && *ldc < f2cmax(1,*n)) { 23175 *info = -13; 23176 } 23177 if (*info != 0) { 23178 i__1 = -(*info); 23179 xerbla_("SBDSQR", &i__1); 23180 return 0; 23181 } 23182 if (*n == 0) { 23183 return 0; 23184 } 23185 if (*n == 1) { 23186 goto L160; 23187 } 23188 23189 /* ROTATE is true if any singular vectors desired, false otherwise */ 23190 23191 rotate = *ncvt > 0 || *nru > 0 || *ncc > 0; 23192 23193 /* If no singular vectors desired, use qd algorithm */ 23194 23195 if (! rotate) { 23196 slasq1_(n, &d__[1], &e[1], &work[1], info); 23197 return 0; 23198 } 23199 23200 nm1 = *n - 1; 23201 nm12 = nm1 + nm1; 23202 nm13 = nm12 + nm1; 23203 idir = 0; 23204 23205 /* Get machine constants */ 23206 23207 eps = slamch_("Epsilon"); 23208 unfl = slamch_("Safe minimum"); 23209 23210 /* If matrix lower bidiagonal, rotate to be upper bidiagonal 23211 by applying Givens rotations on the left */ 23212 23213 if (lower) { 23214 i__1 = *n - 1; 23215 for (i__ = 1; i__ <= i__1; ++i__) { 23216 slartg_(&d__[i__], &e[i__], &cs, &sn, &r__); 23217 d__[i__] = r__; 23218 e[i__] = sn * d__[i__ + 1]; 23219 d__[i__ + 1] = cs * d__[i__ + 1]; 23220 work[i__] = cs; 23221 work[nm1 + i__] = sn; 23222 /* L10: */ 23223 } 23224 23225 /* Update singular vectors if desired */ 23226 23227 if (*nru > 0) { 23228 slasr_("R", "V", "F", nru, n, &work[1], &work[*n], &u[u_offset], 23229 ldu); 23230 } 23231 if (*ncc > 0) { 23232 slasr_("L", "V", "F", n, ncc, &work[1], &work[*n], &c__[c_offset], 23233 ldc); 23234 } 23235 } 23236 23237 /* Compute singular values to relative accuracy TOL 23238 (By setting TOL to be negative, algorithm will compute 23239 singular values to absolute accuracy ABS(TOL)*norm(input matrix)) 23240 23241 Computing MAX 23242 Computing MIN */ 23243 d__1 = (doublereal) eps; 23244 //chao changed pow_dd to pow 23245 r__3 = 100.f, r__4 = pow(d__1, c_b15); 23246 r__1 = 10.f, r__2 = df2cmin(r__3,r__4); 23247 tolmul = df2cmax(r__1,r__2); 23248 tol = tolmul * eps; 23249 23250 /* Compute approximate maximum, minimum singular values */ 23251 23252 smax = 0.f; 23253 i__1 = *n; 23254 for (i__ = 1; i__ <= i__1; ++i__) { 23255 /* Computing MAX */ 23256 r__2 = smax, r__3 = (r__1 = d__[i__], dabs(r__1)); 23257 smax = df2cmax(r__2,r__3); 23258 /* L20: */ 23259 } 23260 i__1 = *n - 1; 23261 for (i__ = 1; i__ <= i__1; ++i__) { 23262 /* Computing MAX */ 23263 r__2 = smax, r__3 = (r__1 = e[i__], dabs(r__1)); 23264 smax = df2cmax(r__2,r__3); 23265 /* L30: */ 23266 } 23267 sminl = 0.f; 23268 if (tol >= 0.f) { 23269 23270 /* Relative accuracy desired */ 23271 23272 sminoa = dabs(d__[1]); 23273 if (sminoa == 0.f) { 23274 goto L50; 23275 } 23276 mu = sminoa; 23277 i__1 = *n; 23278 for (i__ = 2; i__ <= i__1; ++i__) { 23279 mu = (r__2 = d__[i__], dabs(r__2)) * (mu / (mu + (r__1 = e[i__ - 23280 1], dabs(r__1)))); 23281 sminoa = df2cmin(sminoa,mu); 23282 if (sminoa == 0.f) { 23283 goto L50; 23284 } 23285 /* L40: */ 23286 } 23287 L50: 23288 sminoa /= sqrt((real) (*n)); 23289 /* Computing MAX */ 23290 r__1 = tol * sminoa, r__2 = *n * 6 * *n * unfl; 23291 thresh = df2cmax(r__1,r__2); 23292 } else { 23293 23294 /* Absolute accuracy desired 23295 23296 Computing MAX */ 23297 r__1 = dabs(tol) * smax, r__2 = *n * 6 * *n * unfl; 23298 thresh = df2cmax(r__1,r__2); 23299 } 23300 23301 /* Prepare for main iteration loop for the singular values 23302 (MAXIT is the maximum number of passes through the inner 23303 loop permitted before nonconvergence signalled.) */ 23304 23305 maxit = *n * 6 * *n; 23306 iter = 0; 23307 oldll = -1; 23308 oldm = -1; 23309 23310 /* M points to last element of unconverged part of matrix */ 23311 23312 m = *n; 23313 23314 /* Begin main iteration loop */ 23315 23316 L60: 23317 23318 /* Check for convergence or exceeding iteration count */ 23319 23320 if (m <= 1) { 23321 goto L160; 23322 } 23323 if (iter > maxit) { 23324 goto L200; 23325 } 23326 23327 /* Find diagonal block of matrix to work on */ 23328 23329 if (tol < 0.f && (r__1 = d__[m], dabs(r__1)) <= thresh) { 23330 d__[m] = 0.f; 23331 } 23332 smax = (r__1 = d__[m], dabs(r__1)); 23333 smin = smax; 23334 i__1 = m - 1; 23335 for (lll = 1; lll <= i__1; ++lll) { 23336 ll = m - lll; 23337 abss = (r__1 = d__[ll], dabs(r__1)); 23338 abse = (r__1 = e[ll], dabs(r__1)); 23339 if (tol < 0.f && abss <= thresh) { 23340 d__[ll] = 0.f; 23341 } 23342 if (abse <= thresh) { 23343 goto L80; 23344 } 23345 smin = df2cmin(smin,abss); 23346 /* Computing MAX */ 23347 r__1 = f2cmax(smax,abss); 23348 smax = df2cmax(r__1,abse); 23349 /* L70: */ 23350 } 23351 ll = 0; 23352 goto L90; 23353 L80: 23354 e[ll] = 0.f; 23355 23356 /* Matrix splits since E(LL) = 0 */ 23357 23358 if (ll == m - 1) { 23359 23360 /* Convergence of bottom singular value, return to top of loop */ 23361 23362 --m; 23363 goto L60; 23364 } 23365 L90: 23366 ++ll; 23367 23368 /* E(LL) through E(M-1) are nonzero, E(LL-1) is zero */ 23369 23370 if (ll == m - 1) { 23371 23372 /* 2 by 2 block, handle separately */ 23373 23374 slasv2_(&d__[m - 1], &e[m - 1], &d__[m], &sigmn, &sigmx, &sinr, &cosr, 23375 &sinl, &cosl); 23376 d__[m - 1] = sigmx; 23377 e[m - 1] = 0.f; 23378 d__[m] = sigmn; 23379 23380 /* Compute singular vectors, if desired */ 23381 23382 if (*ncvt > 0) { 23383 srot_(ncvt, &vt_ref(m - 1, 1), ldvt, &vt_ref(m, 1), ldvt, &cosr, & 23384 sinr); 23385 } 23386 if (*nru > 0) { 23387 srot_(nru, &u_ref(1, m - 1), &c__1, &u_ref(1, m), &c__1, &cosl, & 23388 sinl); 23389 } 23390 if (*ncc > 0) { 23391 srot_(ncc, &c___ref(m - 1, 1), ldc, &c___ref(m, 1), ldc, &cosl, & 23392 sinl); 23393 } 23394 m += -2; 23395 goto L60; 23396 } 23397 23398 /* If working on new submatrix, choose shift direction 23399 (from larger end diagonal element towards smaller) */ 23400 23401 if (ll > oldm || m < oldll) { 23402 if ((r__1 = d__[ll], dabs(r__1)) >= (r__2 = d__[m], dabs(r__2))) { 23403 23404 /* Chase bulge from top (big end) to bottom (small end) */ 23405 23406 idir = 1; 23407 } else { 23408 23409 /* Chase bulge from bottom (big end) to top (small end) */ 23410 23411 idir = 2; 23412 } 23413 } 23414 23415 /* Apply convergence tests */ 23416 23417 if (idir == 1) { 23418 23419 /* Run convergence test in forward direction 23420 First apply standard test to bottom of matrix */ 23421 23422 if ((r__2 = e[m - 1], dabs(r__2)) <= dabs(tol) * (r__1 = d__[m], dabs( 23423 r__1)) || tol < 0.f && (r__3 = e[m - 1], dabs(r__3)) <= 23424 thresh) { 23425 e[m - 1] = 0.f; 23426 goto L60; 23427 } 23428 23429 if (tol >= 0.f) { 23430 23431 /* If relative accuracy desired, 23432 apply convergence criterion forward */ 23433 23434 mu = (r__1 = d__[ll], dabs(r__1)); 23435 sminl = mu; 23436 i__1 = m - 1; 23437 for (lll = ll; lll <= i__1; ++lll) { 23438 if ((r__1 = e[lll], dabs(r__1)) <= tol * mu) { 23439 e[lll] = 0.f; 23440 goto L60; 23441 } 23442 sminlo = sminl; 23443 mu = (r__2 = d__[lll + 1], dabs(r__2)) * (mu / (mu + (r__1 = 23444 e[lll], dabs(r__1)))); 23445 sminl = df2cmin(sminl,mu); 23446 /* L100: */ 23447 } 23448 } 23449 23450 } else { 23451 23452 /* Run convergence test in backward direction 23453 First apply standard test to top of matrix */ 23454 23455 if ((r__2 = e[ll], dabs(r__2)) <= dabs(tol) * (r__1 = d__[ll], dabs( 23456 r__1)) || tol < 0.f && (r__3 = e[ll], dabs(r__3)) <= thresh) { 23457 e[ll] = 0.f; 23458 goto L60; 23459 } 23460 23461 if (tol >= 0.f) { 23462 23463 /* If relative accuracy desired, 23464 apply convergence criterion backward */ 23465 23466 mu = (r__1 = d__[m], dabs(r__1)); 23467 sminl = mu; 23468 i__1 = ll; 23469 for (lll = m - 1; lll >= i__1; --lll) { 23470 if ((r__1 = e[lll], dabs(r__1)) <= tol * mu) { 23471 e[lll] = 0.f; 23472 goto L60; 23473 } 23474 sminlo = sminl; 23475 mu = (r__2 = d__[lll], dabs(r__2)) * (mu / (mu + (r__1 = e[ 23476 lll], dabs(r__1)))); 23477 sminl = df2cmin(sminl,mu); 23478 /* L110: */ 23479 } 23480 } 23481 } 23482 oldll = ll; 23483 oldm = m; 23484 23485 /* Compute shift. First, test if shifting would ruin relative 23486 accuracy, and if so set the shift to zero. 23487 23488 Computing MAX */ 23489 r__1 = eps, r__2 = tol * .01f; 23490 if (tol >= 0.f && *n * tol * (sminl / smax) <= df2cmax(r__1,r__2)) { 23491 23492 /* Use a zero shift to avoid loss of relative accuracy */ 23493 23494 shift = 0.f; 23495 } else { 23496 23497 /* Compute the shift from 2-by-2 block at end of matrix */ 23498 23499 if (idir == 1) { 23500 sll = (r__1 = d__[ll], dabs(r__1)); 23501 slas2_(&d__[m - 1], &e[m - 1], &d__[m], &shift, &r__); 23502 } else { 23503 sll = (r__1 = d__[m], dabs(r__1)); 23504 slas2_(&d__[ll], &e[ll], &d__[ll + 1], &shift, &r__); 23505 } 23506 23507 /* Test if shift negligible, and if so set to zero */ 23508 23509 if (sll > 0.f) { 23510 /* Computing 2nd power */ 23511 r__1 = shift / sll; 23512 if (r__1 * r__1 < eps) { 23513 shift = 0.f; 23514 } 23515 } 23516 } 23517 23518 /* Increment iteration count */ 23519 23520 iter = iter + m - ll; 23521 23522 /* If SHIFT = 0, do simplified QR iteration */ 23523 23524 if (shift == 0.f) { 23525 if (idir == 1) { 23526 23527 /* Chase bulge from top to bottom 23528 Save cosines and sines for later singular vector updates */ 23529 23530 cs = 1.f; 23531 oldcs = 1.f; 23532 i__1 = m - 1; 23533 for (i__ = ll; i__ <= i__1; ++i__) { 23534 r__1 = d__[i__] * cs; 23535 slartg_(&r__1, &e[i__], &cs, &sn, &r__); 23536 if (i__ > ll) { 23537 e[i__ - 1] = oldsn * r__; 23538 } 23539 r__1 = oldcs * r__; 23540 r__2 = d__[i__ + 1] * sn; 23541 slartg_(&r__1, &r__2, &oldcs, &oldsn, &d__[i__]); 23542 work[i__ - ll + 1] = cs; 23543 work[i__ - ll + 1 + nm1] = sn; 23544 work[i__ - ll + 1 + nm12] = oldcs; 23545 work[i__ - ll + 1 + nm13] = oldsn; 23546 /* L120: */ 23547 } 23548 h__ = d__[m] * cs; 23549 d__[m] = h__ * oldcs; 23550 e[m - 1] = h__ * oldsn; 23551 23552 /* Update singular vectors */ 23553 23554 if (*ncvt > 0) { 23555 i__1 = m - ll + 1; 23556 slasr_("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], & 23557 vt_ref(ll, 1), ldvt); 23558 } 23559 if (*nru > 0) { 23560 i__1 = m - ll + 1; 23561 slasr_("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13 23562 + 1], &u_ref(1, ll), ldu); 23563 } 23564 if (*ncc > 0) { 23565 i__1 = m - ll + 1; 23566 slasr_("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13 23567 + 1], &c___ref(ll, 1), ldc); 23568 } 23569 23570 /* Test convergence */ 23571 23572 if ((r__1 = e[m - 1], dabs(r__1)) <= thresh) { 23573 e[m - 1] = 0.f; 23574 } 23575 23576 } else { 23577 23578 /* Chase bulge from bottom to top 23579 Save cosines and sines for later singular vector updates */ 23580 23581 cs = 1.f; 23582 oldcs = 1.f; 23583 i__1 = ll + 1; 23584 for (i__ = m; i__ >= i__1; --i__) { 23585 r__1 = d__[i__] * cs; 23586 slartg_(&r__1, &e[i__ - 1], &cs, &sn, &r__); 23587 if (i__ < m) { 23588 e[i__] = oldsn * r__; 23589 } 23590 r__1 = oldcs * r__; 23591 r__2 = d__[i__ - 1] * sn; 23592 slartg_(&r__1, &r__2, &oldcs, &oldsn, &d__[i__]); 23593 work[i__ - ll] = cs; 23594 work[i__ - ll + nm1] = -sn; 23595 work[i__ - ll + nm12] = oldcs; 23596 work[i__ - ll + nm13] = -oldsn; 23597 /* L130: */ 23598 } 23599 h__ = d__[ll] * cs; 23600 d__[ll] = h__ * oldcs; 23601 e[ll] = h__ * oldsn; 23602 23603 /* Update singular vectors */ 23604 23605 if (*ncvt > 0) { 23606 i__1 = m - ll + 1; 23607 slasr_("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[ 23608 nm13 + 1], &vt_ref(ll, 1), ldvt); 23609 } 23610 if (*nru > 0) { 23611 i__1 = m - ll + 1; 23612 slasr_("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u_ref( 23613 1, ll), ldu); 23614 } 23615 if (*ncc > 0) { 23616 i__1 = m - ll + 1; 23617 slasr_("L", "V", "B", &i__1, ncc, &work[1], &work[*n], & 23618 c___ref(ll, 1), ldc); 23619 } 23620 23621 /* Test convergence */ 23622 23623 if ((r__1 = e[ll], dabs(r__1)) <= thresh) { 23624 e[ll] = 0.f; 23625 } 23626 } 23627 } else { 23628 23629 /* Use nonzero shift */ 23630 23631 if (idir == 1) { 23632 23633 /* Chase bulge from top to bottom 23634 Save cosines and sines for later singular vector updates */ 23635 23636 f = ((r__1 = d__[ll], dabs(r__1)) - shift) * (r_sign(&c_b49, &d__[ 23637 ll]) + shift / d__[ll]); 23638 g = e[ll]; 23639 i__1 = m - 1; 23640 for (i__ = ll; i__ <= i__1; ++i__) { 23641 slartg_(&f, &g, &cosr, &sinr, &r__); 23642 if (i__ > ll) { 23643 e[i__ - 1] = r__; 23644 } 23645 f = cosr * d__[i__] + sinr * e[i__]; 23646 e[i__] = cosr * e[i__] - sinr * d__[i__]; 23647 g = sinr * d__[i__ + 1]; 23648 d__[i__ + 1] = cosr * d__[i__ + 1]; 23649 slartg_(&f, &g, &cosl, &sinl, &r__); 23650 d__[i__] = r__; 23651 f = cosl * e[i__] + sinl * d__[i__ + 1]; 23652 d__[i__ + 1] = cosl * d__[i__ + 1] - sinl * e[i__]; 23653 if (i__ < m - 1) { 23654 g = sinl * e[i__ + 1]; 23655 e[i__ + 1] = cosl * e[i__ + 1]; 23656 } 23657 work[i__ - ll + 1] = cosr; 23658 work[i__ - ll + 1 + nm1] = sinr; 23659 work[i__ - ll + 1 + nm12] = cosl; 23660 work[i__ - ll + 1 + nm13] = sinl; 23661 /* L140: */ 23662 } 23663 e[m - 1] = f; 23664 23665 /* Update singular vectors */ 23666 23667 if (*ncvt > 0) { 23668 i__1 = m - ll + 1; 23669 slasr_("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], & 23670 vt_ref(ll, 1), ldvt); 23671 } 23672 if (*nru > 0) { 23673 i__1 = m - ll + 1; 23674 slasr_("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13 23675 + 1], &u_ref(1, ll), ldu); 23676 } 23677 if (*ncc > 0) { 23678 i__1 = m - ll + 1; 23679 slasr_("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13 23680 + 1], &c___ref(ll, 1), ldc); 23681 } 23682 23683 /* Test convergence */ 23684 23685 if ((r__1 = e[m - 1], dabs(r__1)) <= thresh) { 23686 e[m - 1] = 0.f; 23687 } 23688 23689 } else { 23690 23691 /* Chase bulge from bottom to top 23692 Save cosines and sines for later singular vector updates */ 23693 23694 f = ((r__1 = d__[m], dabs(r__1)) - shift) * (r_sign(&c_b49, &d__[ 23695 m]) + shift / d__[m]); 23696 g = e[m - 1]; 23697 i__1 = ll + 1; 23698 for (i__ = m; i__ >= i__1; --i__) { 23699 slartg_(&f, &g, &cosr, &sinr, &r__); 23700 if (i__ < m) { 23701 e[i__] = r__; 23702 } 23703 f = cosr * d__[i__] + sinr * e[i__ - 1]; 23704 e[i__ - 1] = cosr * e[i__ - 1] - sinr * d__[i__]; 23705 g = sinr * d__[i__ - 1]; 23706 d__[i__ - 1] = cosr * d__[i__ - 1]; 23707 slartg_(&f, &g, &cosl, &sinl, &r__); 23708 d__[i__] = r__; 23709 f = cosl * e[i__ - 1] + sinl * d__[i__ - 1]; 23710 d__[i__ - 1] = cosl * d__[i__ - 1] - sinl * e[i__ - 1]; 23711 if (i__ > ll + 1) { 23712 g = sinl * e[i__ - 2]; 23713 e[i__ - 2] = cosl * e[i__ - 2]; 23714 } 23715 work[i__ - ll] = cosr; 23716 work[i__ - ll + nm1] = -sinr; 23717 work[i__ - ll + nm12] = cosl; 23718 work[i__ - ll + nm13] = -sinl; 23719 /* L150: */ 23720 } 23721 e[ll] = f; 23722 23723 /* Test convergence */ 23724 23725 if ((r__1 = e[ll], dabs(r__1)) <= thresh) { 23726 e[ll] = 0.f; 23727 } 23728 23729 /* Update singular vectors if desired */ 23730 23731 if (*ncvt > 0) { 23732 i__1 = m - ll + 1; 23733 slasr_("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[ 23734 nm13 + 1], &vt_ref(ll, 1), ldvt); 23735 } 23736 if (*nru > 0) { 23737 i__1 = m - ll + 1; 23738 slasr_("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u_ref( 23739 1, ll), ldu); 23740 } 23741 if (*ncc > 0) { 23742 i__1 = m - ll + 1; 23743 slasr_("L", "V", "B", &i__1, ncc, &work[1], &work[*n], & 23744 c___ref(ll, 1), ldc); 23745 } 23746 } 23747 } 23748 23749 /* QR iteration finished, go back and check convergence */ 23750 23751 goto L60; 23752 23753 /* All singular values converged, so make them positive */ 23754 23755 L160: 23756 i__1 = *n; 23757 for (i__ = 1; i__ <= i__1; ++i__) { 23758 if (d__[i__] < 0.f) { 23759 d__[i__] = -d__[i__]; 23760 23761 /* Change sign of singular vectors, if desired */ 23762 23763 if (*ncvt > 0) { 23764 sscal_(ncvt, &c_b72, &vt_ref(i__, 1), ldvt); 23765 } 23766 } 23767 /* L170: */ 23768 } 23769 23770 /* Sort the singular values into decreasing order (insertion sort on 23771 singular values, but only one transposition per singular vector) */ 23772 23773 i__1 = *n - 1; 23774 for (i__ = 1; i__ <= i__1; ++i__) { 23775 23776 /* Scan for smallest D(I) */ 23777 23778 isub = 1; 23779 smin = d__[1]; 23780 i__2 = *n + 1 - i__; 23781 for (j = 2; j <= i__2; ++j) { 23782 if (d__[j] <= smin) { 23783 isub = j; 23784 smin = d__[j]; 23785 } 23786 /* L180: */ 23787 } 23788 if (isub != *n + 1 - i__) { 23789 23790 /* Swap singular values and vectors */ 23791 23792 d__[isub] = d__[*n + 1 - i__]; 23793 d__[*n + 1 - i__] = smin; 23794 if (*ncvt > 0) { 23795 sswap_(ncvt, &vt_ref(isub, 1), ldvt, &vt_ref(*n + 1 - i__, 1), 23796 ldvt); 23797 } 23798 if (*nru > 0) { 23799 sswap_(nru, &u_ref(1, isub), &c__1, &u_ref(1, *n + 1 - i__), & 23800 c__1); 23801 } 23802 if (*ncc > 0) { 23803 sswap_(ncc, &c___ref(isub, 1), ldc, &c___ref(*n + 1 - i__, 1), 23804 ldc); 23805 } 23806 } 23807 /* L190: */ 23808 } 23809 goto L220; 23810 23811 /* Maximum number of iterations exceeded, failure to converge */ 23812 23813 L200: 23814 *info = 0; 23815 i__1 = *n - 1; 23816 for (i__ = 1; i__ <= i__1; ++i__) { 23817 if (e[i__] != 0.f) { 23818 ++(*info); 23819 } 23820 /* L210: */ 23821 } 23822 L220: 23823 return 0; 23824 23825 /* End of SBDSQR */ 23826 23827 } /* sbdsqr_ */
Definition at line 1163 of file lapackblas.cpp.
References integer.
Referenced by EMAN::PCA::Lanczos_ooc(), slaed0_(), slaed1_(), slaed2_(), slaed3_(), slaed8_(), slaed9_(), slaeda_(), slarfb_(), and slasq1_().
01165 { 01166 /* System generated locals */ 01167 integer i__1; 01168 /* Local variables */ 01169 static integer i__, m, ix, iy, mp1; 01170 /* copies a vector, x, to a vector, y. 01171 uses unrolled loops for increments equal to 1. 01172 jack dongarra, linpack, 3/11/78. 01173 modified 12/3/93, array(1) declarations changed to array(*) 01174 Parameter adjustments */ 01175 --sy; 01176 --sx; 01177 /* Function Body */ 01178 if (*n <= 0) { 01179 return 0; 01180 } 01181 if (*incx == 1 && *incy == 1) { 01182 goto L20; 01183 } 01184 /* code for unequal increments or equal increments 01185 not equal to 1 */ 01186 ix = 1; 01187 iy = 1; 01188 if (*incx < 0) { 01189 ix = (-(*n) + 1) * *incx + 1; 01190 } 01191 if (*incy < 0) { 01192 iy = (-(*n) + 1) * *incy + 1; 01193 } 01194 i__1 = *n; 01195 for (i__ = 1; i__ <= i__1; ++i__) { 01196 sy[iy] = sx[ix]; 01197 ix += *incx; 01198 iy += *incy; 01199 /* L10: */ 01200 } 01201 return 0; 01202 /* code for both increments equal to 1 01203 clean-up loop */ 01204 L20: 01205 m = *n % 7; 01206 if (m == 0) { 01207 goto L40; 01208 } 01209 i__1 = m; 01210 for (i__ = 1; i__ <= i__1; ++i__) { 01211 sy[i__] = sx[i__]; 01212 /* L30: */ 01213 } 01214 if (*n < 7) { 01215 return 0; 01216 } 01217 L40: 01218 mp1 = m + 1; 01219 i__1 = *n; 01220 for (i__ = mp1; i__ <= i__1; i__ += 7) { 01221 sy[i__] = sx[i__]; 01222 sy[i__ + 1] = sx[i__ + 1]; 01223 sy[i__ + 2] = sx[i__ + 2]; 01224 sy[i__ + 3] = sx[i__ + 3]; 01225 sy[i__ + 4] = sx[i__ + 4]; 01226 sy[i__ + 5] = sx[i__ + 5]; 01227 sy[i__ + 6] = sx[i__ + 6]; 01228 /* L50: */ 01229 } 01230 return 0; 01231 } /* scopy_ */
| doublereal sdot_ | ( | integer * | n, | |
| real * | sx, | |||
| integer * | incx, | |||
| real * | sy, | |||
| integer * | incy | |||
| ) |
Definition at line 1236 of file lapackblas.cpp.
References integer.
Referenced by EMAN::PCA::Lanczos(), EMAN::PCAlarge::Lanczos(), EMAN::PCA::Lanczos_ooc(), slatrd_(), and ssytd2_().
01237 { 01238 /* System generated locals */ 01239 integer i__1; 01240 real ret_val; 01241 /* Local variables */ 01242 static integer i__, m; 01243 static real stemp; 01244 static integer ix, iy, mp1; 01245 /* forms the dot product of two vectors. 01246 uses unrolled loops for increments equal to one. 01247 jack dongarra, linpack, 3/11/78. 01248 modified 12/3/93, array(1) declarations changed to array(*) 01249 Parameter adjustments */ 01250 --sy; 01251 --sx; 01252 /* Function Body */ 01253 stemp = 0.f; 01254 ret_val = 0.f; 01255 if (*n <= 0) { 01256 return ret_val; 01257 } 01258 if (*incx == 1 && *incy == 1) { 01259 goto L20; 01260 } 01261 /* code for unequal increments or equal increments 01262 not equal to 1 */ 01263 ix = 1; 01264 iy = 1; 01265 if (*incx < 0) { 01266 ix = (-(*n) + 1) * *incx + 1; 01267 } 01268 if (*incy < 0) { 01269 iy = (-(*n) + 1) * *incy + 1; 01270 } 01271 i__1 = *n; 01272 for (i__ = 1; i__ <= i__1; ++i__) { 01273 stemp += sx[ix] * sy[iy]; 01274 ix += *incx; 01275 iy += *incy; 01276 /* L10: */ 01277 } 01278 ret_val = stemp; 01279 return ret_val; 01280 /* code for both increments equal to 1 01281 clean-up loop */ 01282 L20: 01283 m = *n % 5; 01284 if (m == 0) { 01285 goto L40; 01286 } 01287 i__1 = m; 01288 for (i__ = 1; i__ <= i__1; ++i__) { 01289 stemp += sx[i__] * sy[i__]; 01290 /* L30: */ 01291 } 01292 if (*n < 5) { 01293 goto L60; 01294 } 01295 L40: 01296 mp1 = m + 1; 01297 i__1 = *n; 01298 for (i__ = mp1; i__ <= i__1; i__ += 5) { 01299 stemp = stemp + sx[i__] * sy[i__] + sx[i__ + 1] * sy[i__ + 1] + sx[ 01300 i__ + 2] * sy[i__ + 2] + sx[i__ + 3] * sy[i__ + 3] + sx[i__ + 01301 4] * sy[i__ + 4]; 01302 /* L50: */ 01303 } 01304 L60: 01305 ret_val = stemp; 01306 return ret_val; 01307 } /* sdot_ */
| int sgebd2_ | ( | integer * | m, | |
| integer * | n, | |||
| real * | a, | |||
| integer * | lda, | |||
| real * | d__, | |||
| real * | e, | |||
| real * | tauq, | |||
| real * | taup, | |||
| real * | work, | |||
| integer * | info | |||
| ) |
Definition at line 21343 of file lapackblas.cpp.
References a_ref, f2cmax, f2cmin, integer, slarf_(), slarfg_(), and xerbla_().
Referenced by sgebrd_().
21345 { 21346 /* -- LAPACK routine (version 3.0) -- 21347 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 21348 Courant Institute, Argonne National Lab, and Rice University 21349 February 29, 1992 21350 21351 21352 Purpose 21353 ======= 21354 21355 SGEBD2 reduces a real general m by n matrix A to upper or lower 21356 bidiagonal form B by an orthogonal transformation: Q' * A * P = B. 21357 21358 If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. 21359 21360 Arguments 21361 ========= 21362 21363 M (input) INTEGER 21364 The number of rows in the matrix A. M >= 0. 21365 21366 N (input) INTEGER 21367 The number of columns in the matrix A. N >= 0. 21368 21369 A (input/output) REAL array, dimension (LDA,N) 21370 On entry, the m by n general matrix to be reduced. 21371 On exit, 21372 if m >= n, the diagonal and the first superdiagonal are 21373 overwritten with the upper bidiagonal matrix B; the 21374 elements below the diagonal, with the array TAUQ, represent 21375 the orthogonal matrix Q as a product of elementary 21376 reflectors, and the elements above the first superdiagonal, 21377 with the array TAUP, represent the orthogonal matrix P as 21378 a product of elementary reflectors; 21379 if m < n, the diagonal and the first subdiagonal are 21380 overwritten with the lower bidiagonal matrix B; the 21381 elements below the first subdiagonal, with the array TAUQ, 21382 represent the orthogonal matrix Q as a product of 21383 elementary reflectors, and the elements above the diagonal, 21384 with the array TAUP, represent the orthogonal matrix P as 21385 a product of elementary reflectors. 21386 See Further Details. 21387 21388 LDA (input) INTEGER 21389 The leading dimension of the array A. LDA >= max(1,M). 21390 21391 D (output) REAL array, dimension (min(M,N)) 21392 The diagonal elements of the bidiagonal matrix B: 21393 D(i) = A(i,i). 21394 21395 E (output) REAL array, dimension (min(M,N)-1) 21396 The off-diagonal elements of the bidiagonal matrix B: 21397 if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1; 21398 if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1. 21399 21400 TAUQ (output) REAL array dimension (min(M,N)) 21401 The scalar factors of the elementary reflectors which 21402 represent the orthogonal matrix Q. See Further Details. 21403 21404 TAUP (output) REAL array, dimension (min(M,N)) 21405 The scalar factors of the elementary reflectors which 21406 represent the orthogonal matrix P. See Further Details. 21407 21408 WORK (workspace) REAL array, dimension (max(M,N)) 21409 21410 INFO (output) INTEGER 21411 = 0: successful exit. 21412 < 0: if INFO = -i, the i-th argument had an illegal value. 21413 21414 Further Details 21415 =============== 21416 21417 The matrices Q and P are represented as products of elementary 21418 reflectors: 21419 21420 If m >= n, 21421 21422 Q = H(1) H(2) . . . H(n) and P = G(1) G(2) . . . G(n-1) 21423 21424 Each H(i) and G(i) has the form: 21425 21426 H(i) = I - tauq * v * v' and G(i) = I - taup * u * u' 21427 21428 where tauq and taup are real scalars, and v and u are real vectors; 21429 v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i); 21430 u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n); 21431 tauq is stored in TAUQ(i) and taup in TAUP(i). 21432 21433 If m < n, 21434 21435 Q = H(1) H(2) . . . H(m-1) and P = G(1) G(2) . . . G(m) 21436 21437 Each H(i) and G(i) has the form: 21438 21439 H(i) = I - tauq * v * v' and G(i) = I - taup * u * u' 21440 21441 where tauq and taup are real scalars, and v and u are real vectors; 21442 v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i); 21443 u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n); 21444 tauq is stored in TAUQ(i) and taup in TAUP(i). 21445 21446 The contents of A on exit are illustrated by the following examples: 21447 21448 m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n): 21449 21450 ( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 ) 21451 ( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 ) 21452 ( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 ) 21453 ( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 ) 21454 ( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 ) 21455 ( v1 v2 v3 v4 v5 ) 21456 21457 where d and e denote diagonal and off-diagonal elements of B, vi 21458 denotes an element of the vector defining H(i), and ui an element of 21459 the vector defining G(i). 21460 21461 ===================================================================== 21462 21463 21464 Test the input parameters 21465 21466 Parameter adjustments */ 21467 /* Table of constant values */ 21468 static integer c__1 = 1; 21469 21470 /* System generated locals */ 21471 integer a_dim1, a_offset, i__1, i__2, i__3, i__4; 21472 /* Local variables */ 21473 static integer i__; 21474 extern /* Subroutine */ int slarf_(const char *, integer *, integer *, real *, 21475 integer *, real *, real *, integer *, real *), xerbla_( 21476 const char *, integer *), slarfg_(integer *, real *, real *, 21477 integer *, real *); 21478 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 21479 21480 21481 a_dim1 = *lda; 21482 a_offset = 1 + a_dim1 * 1; 21483 a -= a_offset; 21484 --d__; 21485 --e; 21486 --tauq; 21487 --taup; 21488 --work; 21489 21490 /* Function Body */ 21491 *info = 0; 21492 if (*m < 0) { 21493 *info = -1; 21494 } else if (*n < 0) { 21495 *info = -2; 21496 } else if (*lda < f2cmax(1,*m)) { 21497 *info = -4; 21498 } 21499 if (*info < 0) { 21500 i__1 = -(*info); 21501 xerbla_("SGEBD2", &i__1); 21502 return 0; 21503 } 21504 21505 if (*m >= *n) { 21506 21507 /* Reduce to upper bidiagonal form */ 21508 21509 i__1 = *n; 21510 for (i__ = 1; i__ <= i__1; ++i__) { 21511 21512 /* Generate elementary reflector H(i) to annihilate A(i+1:m,i) 21513 21514 Computing MIN */ 21515 i__2 = i__ + 1; 21516 i__3 = *m - i__ + 1; 21517 slarfg_(&i__3, &a_ref(i__, i__), &a_ref(f2cmin(i__2,*m), i__), &c__1, 21518 &tauq[i__]); 21519 d__[i__] = a_ref(i__, i__); 21520 a_ref(i__, i__) = 1.f; 21521 21522 /* Apply H(i) to A(i:m,i+1:n) from the left */ 21523 21524 i__2 = *m - i__ + 1; 21525 i__3 = *n - i__; 21526 slarf_("Left", &i__2, &i__3, &a_ref(i__, i__), &c__1, &tauq[i__], 21527 &a_ref(i__, i__ + 1), lda, &work[1]); 21528 a_ref(i__, i__) = d__[i__]; 21529 21530 if (i__ < *n) { 21531 21532 /* Generate elementary reflector G(i) to annihilate 21533 A(i,i+2:n) 21534 21535 Computing MIN */ 21536 i__2 = i__ + 2; 21537 i__3 = *n - i__; 21538 slarfg_(&i__3, &a_ref(i__, i__ + 1), &a_ref(i__, f2cmin(i__2,*n)) 21539 , lda, &taup[i__]); 21540 e[i__] = a_ref(i__, i__ + 1); 21541 a_ref(i__, i__ + 1) = 1.f; 21542 21543 /* Apply G(i) to A(i+1:m,i+1:n) from the right */ 21544 21545 i__2 = *m - i__; 21546 i__3 = *n - i__; 21547 slarf_("Right", &i__2, &i__3, &a_ref(i__, i__ + 1), lda, & 21548 taup[i__], &a_ref(i__ + 1, i__ + 1), lda, &work[1]); 21549 a_ref(i__, i__ + 1) = e[i__]; 21550 } else { 21551 taup[i__] = 0.f; 21552 } 21553 /* L10: */ 21554 } 21555 } else { 21556 21557 /* Reduce to lower bidiagonal form */ 21558 21559 i__1 = *m; 21560 for (i__ = 1; i__ <= i__1; ++i__) { 21561 21562 /* Generate elementary reflector G(i) to annihilate A(i,i+1:n) 21563 21564 Computing MIN */ 21565 i__2 = i__ + 1; 21566 i__3 = *n - i__ + 1; 21567 slarfg_(&i__3, &a_ref(i__, i__), &a_ref(i__, f2cmin(i__2,*n)), lda, & 21568 taup[i__]); 21569 d__[i__] = a_ref(i__, i__); 21570 a_ref(i__, i__) = 1.f; 21571 21572 /* Apply G(i) to A(i+1:m,i:n) from the right 21573 21574 Computing MIN */ 21575 i__2 = i__ + 1; 21576 i__3 = *m - i__; 21577 i__4 = *n - i__ + 1; 21578 slarf_("Right", &i__3, &i__4, &a_ref(i__, i__), lda, &taup[i__], & 21579 a_ref(f2cmin(i__2,*m), i__), lda, &work[1]); 21580 a_ref(i__, i__) = d__[i__]; 21581 21582 if (i__ < *m) { 21583 21584 /* Generate elementary reflector H(i) to annihilate 21585 A(i+2:m,i) 21586 21587 Computing MIN */ 21588 i__2 = i__ + 2; 21589 i__3 = *m - i__; 21590 slarfg_(&i__3, &a_ref(i__ + 1, i__), &a_ref(f2cmin(i__2,*m), i__) 21591 , &c__1, &tauq[i__]); 21592 e[i__] = a_ref(i__ + 1, i__); 21593 a_ref(i__ + 1, i__) = 1.f; 21594 21595 /* Apply H(i) to A(i+1:m,i+1:n) from the left */ 21596 21597 i__2 = *m - i__; 21598 i__3 = *n - i__; 21599 slarf_("Left", &i__2, &i__3, &a_ref(i__ + 1, i__), &c__1, & 21600 tauq[i__], &a_ref(i__ + 1, i__ + 1), lda, &work[1]); 21601 a_ref(i__ + 1, i__) = e[i__]; 21602 } else { 21603 tauq[i__] = 0.f; 21604 } 21605 /* L20: */ 21606 } 21607 } 21608 return 0; 21609 21610 /* End of SGEBD2 */ 21611 21612 } /* sgebd2_ */
| int sgebrd_ | ( | integer * | m, | |
| integer * | n, | |||
| real * | a, | |||
| integer * | lda, | |||
| real * | d__, | |||
| real * | e, | |||
| real * | tauq, | |||
| real * | taup, | |||
| real * | work, | |||
| integer * | lwork, | |||
| integer * | info | |||
| ) |
Definition at line 21027 of file lapackblas.cpp.
References a_ref, c__3, f2cmax, f2cmin, ilaenv_(), integer, nx, sgebd2_(), sgemm_(), slabrd_(), and xerbla_().
Referenced by sgesvd_().
21030 { 21031 /* -- LAPACK routine (version 3.0) -- 21032 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 21033 Courant Institute, Argonne National Lab, and Rice University 21034 June 30, 1999 21035 21036 21037 Purpose 21038 ======= 21039 21040 SGEBRD reduces a general real M-by-N matrix A to upper or lower 21041 bidiagonal form B by an orthogonal transformation: Q**T * A * P = B. 21042 21043 If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. 21044 21045 Arguments 21046 ========= 21047 21048 M (input) INTEGER 21049 The number of rows in the matrix A. M >= 0. 21050 21051 N (input) INTEGER 21052 The number of columns in the matrix A. N >= 0. 21053 21054 A (input/output) REAL array, dimension (LDA,N) 21055 On entry, the M-by-N general matrix to be reduced. 21056 On exit, 21057 if m >= n, the diagonal and the first superdiagonal are 21058 overwritten with the upper bidiagonal matrix B; the 21059 elements below the diagonal, with the array TAUQ, represent 21060 the orthogonal matrix Q as a product of elementary 21061 reflectors, and the elements above the first superdiagonal, 21062 with the array TAUP, represent the orthogonal matrix P as 21063 a product of elementary reflectors; 21064 if m < n, the diagonal and the first subdiagonal are 21065 overwritten with the lower bidiagonal matrix B; the 21066 elements below the first subdiagonal, with the array TAUQ, 21067 represent the orthogonal matrix Q as a product of 21068 elementary reflectors, and the elements above the diagonal, 21069 with the array TAUP, represent the orthogonal matrix P as 21070 a product of elementary reflectors. 21071 See Further Details. 21072 21073 LDA (input) INTEGER 21074 The leading dimension of the array A. LDA >= max(1,M). 21075 21076 D (output) REAL array, dimension (min(M,N)) 21077 The diagonal elements of the bidiagonal matrix B: 21078 D(i) = A(i,i). 21079 21080 E (output) REAL array, dimension (min(M,N)-1) 21081 The off-diagonal elements of the bidiagonal matrix B: 21082 if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1; 21083 if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1. 21084 21085 TAUQ (output) REAL array dimension (min(M,N)) 21086 The scalar factors of the elementary reflectors which 21087 represent the orthogonal matrix Q. See Further Details. 21088 21089 TAUP (output) REAL array, dimension (min(M,N)) 21090 The scalar factors of the elementary reflectors which 21091 represent the orthogonal matrix P. See Further Details. 21092 21093 WORK (workspace/output) REAL array, dimension (LWORK) 21094 On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 21095 21096 LWORK (input) INTEGER 21097 The length of the array WORK. LWORK >= max(1,M,N). 21098 For optimum performance LWORK >= (M+N)*NB, where NB 21099 is the optimal blocksize. 21100 21101 If LWORK = -1, then a workspace query is assumed; the routine 21102 only calculates the optimal size of the WORK array, returns 21103 this value as the first entry of the WORK array, and no error 21104 message related to LWORK is issued by XERBLA. 21105 21106 INFO (output) INTEGER 21107 = 0: successful exit 21108 < 0: if INFO = -i, the i-th argument had an illegal value. 21109 21110 Further Details 21111 =============== 21112 21113 The matrices Q and P are represented as products of elementary 21114 reflectors: 21115 21116 If m >= n, 21117 21118 Q = H(1) H(2) . . . H(n) and P = G(1) G(2) . . . G(n-1) 21119 21120 Each H(i) and G(i) has the form: 21121 21122 H(i) = I - tauq * v * v' and G(i) = I - taup * u * u' 21123 21124 where tauq and taup are real scalars, and v and u are real vectors; 21125 v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i); 21126 u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n); 21127 tauq is stored in TAUQ(i) and taup in TAUP(i). 21128 21129 If m < n, 21130 21131 Q = H(1) H(2) . . . H(m-1) and P = G(1) G(2) . . . G(m) 21132 21133 Each H(i) and G(i) has the form: 21134 21135 H(i) = I - tauq * v * v' and G(i) = I - taup * u * u' 21136 21137 where tauq and taup are real scalars, and v and u are real vectors; 21138 v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i); 21139 u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n); 21140 tauq is stored in TAUQ(i) and taup in TAUP(i). 21141 21142 The contents of A on exit are illustrated by the following examples: 21143 21144 m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n): 21145 21146 ( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 ) 21147 ( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 ) 21148 ( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 ) 21149 ( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 ) 21150 ( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 ) 21151 ( v1 v2 v3 v4 v5 ) 21152 21153 where d and e denote diagonal and off-diagonal elements of B, vi 21154 denotes an element of the vector defining H(i), and ui an element of 21155 the vector defining G(i). 21156 21157 ===================================================================== 21158 21159 21160 Test the input parameters 21161 21162 Parameter adjustments */ 21163 /* Table of constant values */ 21164 static integer c__1 = 1; 21165 static integer c_n1 = -1; 21166 static integer c__3 = 3; 21167 static integer c__2 = 2; 21168 static real c_b21 = -1.f; 21169 static real c_b22 = 1.f; 21170 21171 /* System generated locals */ 21172 integer a_dim1, a_offset, i__1, i__2, i__3, i__4; 21173 /* Local variables */ 21174 static integer i__, j, nbmin, iinfo; 21175 extern /* Subroutine */ int sgemm_(const char *, const char *, integer *, integer *, 21176 integer *, real *, real *, integer *, real *, integer *, real *, 21177 real *, integer *); 21178 static integer minmn; 21179 extern /* Subroutine */ int sgebd2_(integer *, integer *, real *, integer 21180 *, real *, real *, real *, real *, real *, integer *); 21181 static integer nb, nx; 21182 extern /* Subroutine */ int slabrd_(integer *, integer *, integer *, real 21183 *, integer *, real *, real *, real *, real *, real *, integer *, 21184 real *, integer *); 21185 static real ws; 21186 extern /* Subroutine */ int xerbla_(const char *, integer *); 21187 extern integer ilaenv_(integer *, const char *, const char *, integer *, integer *, 21188 integer *, integer *, ftnlen, ftnlen); 21189 static integer ldwrkx, ldwrky, lwkopt; 21190 static logical lquery; 21191 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 21192 21193 21194 a_dim1 = *lda; 21195 a_offset = 1 + a_dim1 * 1; 21196 a -= a_offset; 21197 --d__; 21198 --e; 21199 --tauq; 21200 --taup; 21201 --work; 21202 21203 /* Function Body */ 21204 *info = 0; 21205 /* Computing MAX */ 21206 i__1 = 1, i__2 = ilaenv_(&c__1, "SGEBRD", " ", m, n, &c_n1, &c_n1, ( 21207 ftnlen)6, (ftnlen)1); 21208 nb = f2cmax(i__1,i__2); 21209 lwkopt = (*m + *n) * nb; 21210 work[1] = (real) lwkopt; 21211 lquery = *lwork == -1; 21212 if (*m < 0) { 21213 *info = -1; 21214 } else if (*n < 0) { 21215 *info = -2; 21216 } else if (*lda < f2cmax(1,*m)) { 21217 *info = -4; 21218 } else /* if(complicated condition) */ { 21219 /* Computing MAX */ 21220 i__1 = f2cmax(1,*m); 21221 if (*lwork < f2cmax(i__1,*n) && ! lquery) { 21222 *info = -10; 21223 } 21224 } 21225 if (*info < 0) { 21226 i__1 = -(*info); 21227 xerbla_("SGEBRD", &i__1); 21228 return 0; 21229 } else if (lquery) { 21230 return 0; 21231 } 21232 21233 /* Quick return if possible */ 21234 21235 minmn = f2cmin(*m,*n); 21236 if (minmn == 0) { 21237 work[1] = 1.f; 21238 return 0; 21239 } 21240 21241 ws = (real) f2cmax(*m,*n); 21242 ldwrkx = *m; 21243 ldwrky = *n; 21244 21245 if (nb > 1 && nb < minmn) { 21246 21247 /* Set the crossover point NX. 21248 21249 Computing MAX */ 21250 i__1 = nb, i__2 = ilaenv_(&c__3, "SGEBRD", " ", m, n, &c_n1, &c_n1, ( 21251 ftnlen)6, (ftnlen)1); 21252 nx = f2cmax(i__1,i__2); 21253 21254 /* Determine when to switch from blocked to unblocked code. */ 21255 21256 if (nx < minmn) { 21257 ws = (real) ((*m + *n) * nb); 21258 if ((real) (*lwork) < ws) { 21259 21260 /* Not enough work space for the optimal NB, consider using 21261 a smaller block size. */ 21262 21263 nbmin = ilaenv_(&c__2, "SGEBRD", " ", m, n, &c_n1, &c_n1, ( 21264 ftnlen)6, (ftnlen)1); 21265 if (*lwork >= (*m + *n) * nbmin) { 21266 nb = *lwork / (*m + *n); 21267 } else { 21268 nb = 1; 21269 nx = minmn; 21270 } 21271 } 21272 } 21273 } else { 21274 nx = minmn; 21275 } 21276 21277 i__1 = minmn - nx; 21278 i__2 = nb; 21279 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) { 21280 21281 /* Reduce rows and columns i:i+nb-1 to bidiagonal form and return 21282 the matrices X and Y which are needed to update the unreduced 21283 part of the matrix */ 21284 21285 i__3 = *m - i__ + 1; 21286 i__4 = *n - i__ + 1; 21287 slabrd_(&i__3, &i__4, &nb, &a_ref(i__, i__), lda, &d__[i__], &e[i__], 21288 &tauq[i__], &taup[i__], &work[1], &ldwrkx, &work[ldwrkx * nb 21289 + 1], &ldwrky); 21290 21291 /* Update the trailing submatrix A(i+nb:m,i+nb:n), using an update 21292 of the form A := A - V*Y' - X*U' */ 21293 21294 i__3 = *m - i__ - nb + 1; 21295 i__4 = *n - i__ - nb + 1; 21296 sgemm_("No transpose", "Transpose", &i__3, &i__4, &nb, &c_b21, &a_ref( 21297 i__ + nb, i__), lda, &work[ldwrkx * nb + nb + 1], &ldwrky, & 21298 c_b22, &a_ref(i__ + nb, i__ + nb), lda) 21299 ; 21300 i__3 = *m - i__ - nb + 1; 21301 i__4 = *n - i__ - nb + 1; 21302 sgemm_("No transpose", "No transpose", &i__3, &i__4, &nb, &c_b21, & 21303 work[nb + 1], &ldwrkx, &a_ref(i__, i__ + nb), lda, &c_b22, & 21304 a_ref(i__ + nb, i__ + nb), lda); 21305 21306 /* Copy diagonal and off-diagonal elements of B back into A */ 21307 21308 if (*m >= *n) { 21309 i__3 = i__ + nb - 1; 21310 for (j = i__; j <= i__3; ++j) { 21311 a_ref(j, j) = d__[j]; 21312 a_ref(j, j + 1) = e[j]; 21313 /* L10: */ 21314 } 21315 } else { 21316 i__3 = i__ + nb - 1; 21317 for (j = i__; j <= i__3; ++j) { 21318 a_ref(j, j) = d__[j]; 21319 a_ref(j + 1, j) = e[j]; 21320 /* L20: */ 21321 } 21322 } 21323 /* L30: */ 21324 } 21325 21326 /* Use unblocked code to reduce the remainder of the matrix */ 21327 21328 i__2 = *m - i__ + 1; 21329 i__1 = *n - i__ + 1; 21330 sgebd2_(&i__2, &i__1, &a_ref(i__, i__), lda, &d__[i__], &e[i__], &tauq[ 21331 i__], &taup[i__], &work[1], &iinfo); 21332 work[1] = ws; 21333 return 0; 21334 21335 /* End of SGEBRD */ 21336 21337 } /* sgebrd_ */
| int sgelq2_ | ( | integer * | m, | |
| integer * | n, | |||
| real * | a, | |||
| integer * | lda, | |||
| real * | tau, | |||
| real * | work, | |||
| integer * | info | |||
| ) |
Definition at line 22808 of file lapackblas.cpp.
References a_ref, f2cmax, f2cmin, integer, slarf_(), slarfg_(), and xerbla_().
Referenced by sgelqf_().
22810 { 22811 /* -- LAPACK routine (version 3.0) -- 22812 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 22813 Courant Institute, Argonne National Lab, and Rice University 22814 February 29, 1992 22815 22816 22817 Purpose 22818 ======= 22819 22820 SGELQ2 computes an LQ factorization of a real m by n matrix A: 22821 A = L * Q. 22822 22823 Arguments 22824 ========= 22825 22826 M (input) INTEGER 22827 The number of rows of the matrix A. M >= 0. 22828 22829 N (input) INTEGER 22830 The number of columns of the matrix A. N >= 0. 22831 22832 A (input/output) REAL array, dimension (LDA,N) 22833 On entry, the m by n matrix A. 22834 On exit, the elements on and below the diagonal of the array 22835 contain the m by min(m,n) lower trapezoidal matrix L (L is 22836 lower triangular if m <= n); the elements above the diagonal, 22837 with the array TAU, represent the orthogonal matrix Q as a 22838 product of elementary reflectors (see Further Details). 22839 22840 LDA (input) INTEGER 22841 The leading dimension of the array A. LDA >= max(1,M). 22842 22843 TAU (output) REAL array, dimension (min(M,N)) 22844 The scalar factors of the elementary reflectors (see Further 22845 Details). 22846 22847 WORK (workspace) REAL array, dimension (M) 22848 22849 INFO (output) INTEGER 22850 = 0: successful exit 22851 < 0: if INFO = -i, the i-th argument had an illegal value 22852 22853 Further Details 22854 =============== 22855 22856 The matrix Q is represented as a product of elementary reflectors 22857 22858 Q = H(k) . . . H(2) H(1), where k = min(m,n). 22859 22860 Each H(i) has the form 22861 22862 H(i) = I - tau * v * v' 22863 22864 where tau is a real scalar, and v is a real vector with 22865 v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i,i+1:n), 22866 and tau in TAU(i). 22867 22868 ===================================================================== 22869 22870 22871 Test the input arguments 22872 22873 Parameter adjustments */ 22874 /* System generated locals */ 22875 integer a_dim1, a_offset, i__1, i__2, i__3; 22876 /* Local variables */ 22877 static integer i__, k; 22878 extern /* Subroutine */ int slarf_(const char *, integer *, integer *, real *, 22879 integer *, real *, real *, integer *, real *), xerbla_( 22880 const char *, integer *), slarfg_(integer *, real *, real *, 22881 integer *, real *); 22882 static real aii; 22883 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 22884 22885 a_dim1 = *lda; 22886 a_offset = 1 + a_dim1 * 1; 22887 a -= a_offset; 22888 --tau; 22889 --work; 22890 22891 /* Function Body */ 22892 *info = 0; 22893 if (*m < 0) { 22894 *info = -1; 22895 } else if (*n < 0) { 22896 *info = -2; 22897 } else if (*lda < f2cmax(1,*m)) { 22898 *info = -4; 22899 } 22900 if (*info != 0) { 22901 i__1 = -(*info); 22902 xerbla_("SGELQ2", &i__1); 22903 return 0; 22904 } 22905 22906 k = f2cmin(*m,*n); 22907 22908 i__1 = k; 22909 for (i__ = 1; i__ <= i__1; ++i__) { 22910 22911 /* Generate elementary reflector H(i) to annihilate A(i,i+1:n) 22912 22913 Computing MIN */ 22914 i__2 = i__ + 1; 22915 i__3 = *n - i__ + 1; 22916 slarfg_(&i__3, &a_ref(i__, i__), &a_ref(i__, f2cmin(i__2,*n)), lda, &tau[ 22917 i__]); 22918 if (i__ < *m) { 22919 22920 /* Apply H(i) to A(i+1:m,i:n) from the right */ 22921 22922 aii = a_ref(i__, i__); 22923 a_ref(i__, i__) = 1.f; 22924 i__2 = *m - i__; 22925 i__3 = *n - i__ + 1; 22926 slarf_("Right", &i__2, &i__3, &a_ref(i__, i__), lda, &tau[i__], & 22927 a_ref(i__ + 1, i__), lda, &work[1]); 22928 a_ref(i__, i__) = aii; 22929 } 22930 /* L10: */ 22931 } 22932 return 0; 22933 22934 /* End of SGELQ2 */ 22935 22936 } /* sgelq2_ */
| int sgelqf_ | ( | integer * | m, | |
| integer * | n, | |||
| real * | a, | |||
| integer * | lda, | |||
| real * | tau, | |||
| real * | work, | |||
| integer * | lwork, | |||
| integer * | info | |||
| ) |
Definition at line 21958 of file lapackblas.cpp.
References a_ref, c__3, f2cmax, f2cmin, ilaenv_(), integer, nx, sgelq2_(), slarfb_(), slarft_(), and xerbla_().
Referenced by sgesvd_().
21960 { 21961 /* -- LAPACK routine (version 3.0) -- 21962 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 21963 Courant Institute, Argonne National Lab, and Rice University 21964 June 30, 1999 21965 21966 21967 Purpose 21968 ======= 21969 21970 SGELQF computes an LQ factorization of a real M-by-N matrix A: 21971 A = L * Q. 21972 21973 Arguments 21974 ========= 21975 21976 M (input) INTEGER 21977 The number of rows of the matrix A. M >= 0. 21978 21979 N (input) INTEGER 21980 The number of columns of the matrix A. N >= 0. 21981 21982 A (input/output) REAL array, dimension (LDA,N) 21983 On entry, the M-by-N matrix A. 21984 On exit, the elements on and below the diagonal of the array 21985 contain the m-by-min(m,n) lower trapezoidal matrix L (L is 21986 lower triangular if m <= n); the elements above the diagonal, 21987 with the array TAU, represent the orthogonal matrix Q as a 21988 product of elementary reflectors (see Further Details). 21989 21990 LDA (input) INTEGER 21991 The leading dimension of the array A. LDA >= max(1,M). 21992 21993 TAU (output) REAL array, dimension (min(M,N)) 21994 The scalar factors of the elementary reflectors (see Further 21995 Details). 21996 21997 WORK (workspace/output) REAL array, dimension (LWORK) 21998 On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 21999 22000 LWORK (input) INTEGER 22001 The dimension of the array WORK. LWORK >= max(1,M). 22002 For optimum performance LWORK >= M*NB, where NB is the 22003 optimal blocksize. 22004 22005 If LWORK = -1, then a workspace query is assumed; the routine 22006 only calculates the optimal size of the WORK array, returns 22007 this value as the first entry of the WORK array, and no error 22008 message related to LWORK is issued by XERBLA. 22009 22010 INFO (output) INTEGER 22011 = 0: successful exit 22012 < 0: if INFO = -i, the i-th argument had an illegal value 22013 22014 Further Details 22015 =============== 22016 22017 The matrix Q is represented as a product of elementary reflectors 22018 22019 Q = H(k) . . . H(2) H(1), where k = min(m,n). 22020 22021 Each H(i) has the form 22022 22023 H(i) = I - tau * v * v' 22024 22025 where tau is a real scalar, and v is a real vector with 22026 v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i,i+1:n), 22027 and tau in TAU(i). 22028 22029 ===================================================================== 22030 22031 22032 Test the input arguments 22033 22034 Parameter adjustments */ 22035 /* Table of constant values */ 22036 static integer c__1 = 1; 22037 static integer c_n1 = -1; 22038 static integer c__3 = 3; 22039 static integer c__2 = 2; 22040 22041 /* System generated locals */ 22042 integer a_dim1, a_offset, i__1, i__2, i__3, i__4; 22043 /* Local variables */ 22044 static integer i__, k, nbmin, iinfo; 22045 extern /* Subroutine */ int sgelq2_(integer *, integer *, real *, integer 22046 *, real *, real *, integer *); 22047 static integer ib, nb, nx; 22048 extern /* Subroutine */ int slarfb_(const char *, const char *, const char *, const char *, 22049 integer *, integer *, integer *, real *, integer *, real *, 22050 integer *, real *, integer *, real *, integer *), xerbla_(const char *, integer *); 22051 extern integer ilaenv_(integer *, const char *, const char *, integer *, integer *, 22052 integer *, integer *, ftnlen, ftnlen); 22053 extern /* Subroutine */ int slarft_(const char *, const char *, integer *, integer *, 22054 real *, integer *, real *, real *, integer *); 22055 static integer ldwork, lwkopt; 22056 static logical lquery; 22057 static integer iws; 22058 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 22059 22060 22061 a_dim1 = *lda; 22062 a_offset = 1 + a_dim1 * 1; 22063 a -= a_offset; 22064 --tau; 22065 --work; 22066 22067 /* Function Body */ 22068 *info = 0; 22069 nb = ilaenv_(&c__1, "SGELQF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen) 22070 1); 22071 lwkopt = *m * nb; 22072 work[1] = (real) lwkopt; 22073 lquery = *lwork == -1; 22074 if (*m < 0) { 22075 *info = -1; 22076 } else if (*n < 0) { 22077 *info = -2; 22078 } else if (*lda < f2cmax(1,*m)) { 22079 *info = -4; 22080 } else if (*lwork < f2cmax(1,*m) && ! lquery) { 22081 *info = -7; 22082 } 22083 if (*info != 0) { 22084 i__1 = -(*info); 22085 xerbla_("SGELQF", &i__1); 22086 return 0; 22087 } else if (lquery) { 22088 return 0; 22089 } 22090 22091 /* Quick return if possible */ 22092 22093 k = f2cmin(*m,*n); 22094 if (k == 0) { 22095 work[1] = 1.f; 22096 return 0; 22097 } 22098 22099 nbmin = 2; 22100 nx = 0; 22101 iws = *m; 22102 if (nb > 1 && nb < k) { 22103 22104 /* Determine when to cross over from blocked to unblocked code. 22105 22106 Computing MAX */ 22107 i__1 = 0, i__2 = ilaenv_(&c__3, "SGELQF", " ", m, n, &c_n1, &c_n1, ( 22108 ftnlen)6, (ftnlen)1); 22109 nx = f2cmax(i__1,i__2); 22110 if (nx < k) { 22111 22112 /* Determine if workspace is large enough for blocked code. */ 22113 22114 ldwork = *m; 22115 iws = ldwork * nb; 22116 if (*lwork < iws) { 22117 22118 /* Not enough workspace to use optimal NB: reduce NB and 22119 determine the minimum value of NB. */ 22120 22121 nb = *lwork / ldwork; 22122 /* Computing MAX */ 22123 i__1 = 2, i__2 = ilaenv_(&c__2, "SGELQF", " ", m, n, &c_n1, & 22124 c_n1, (ftnlen)6, (ftnlen)1); 22125 nbmin = f2cmax(i__1,i__2); 22126 } 22127 } 22128 } 22129 22130 if (nb >= nbmin && nb < k && nx < k) { 22131 22132 /* Use blocked code initially */ 22133 22134 i__1 = k - nx; 22135 i__2 = nb; 22136 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) { 22137 /* Computing MIN */ 22138 i__3 = k - i__ + 1; 22139 ib = f2cmin(i__3,nb); 22140 22141 /* Compute the LQ factorization of the current block 22142 A(i:i+ib-1,i:n) */ 22143 22144 i__3 = *n - i__ + 1; 22145 sgelq2_(&ib, &i__3, &a_ref(i__, i__), lda, &tau[i__], &work[1], & 22146 iinfo); 22147 if (i__ + ib <= *m) { 22148 22149 /* Form the triangular factor of the block reflector 22150 H = H(i) H(i+1) . . . H(i+ib-1) */ 22151 22152 i__3 = *n - i__ + 1; 22153 slarft_("Forward", "Rowwise", &i__3, &ib, &a_ref(i__, i__), 22154 lda, &tau[i__], &work[1], &ldwork); 22155 22156 /* Apply H to A(i+ib:m,i:n) from the right */ 22157 22158 i__3 = *m - i__ - ib + 1; 22159 i__4 = *n - i__ + 1; 22160 slarfb_("Right", "No transpose", "Forward", "Rowwise", &i__3, 22161 &i__4, &ib, &a_ref(i__, i__), lda, &work[1], &ldwork, 22162 &a_ref(i__ + ib, i__), lda, &work[ib + 1], &ldwork); 22163 } 22164 /* L10: */ 22165 } 22166 } else { 22167 i__ = 1; 22168 } 22169 22170 /* Use unblocked code to factor the last or only block. */ 22171 22172 if (i__ <= k) { 22173 i__2 = *m - i__ + 1; 22174 i__1 = *n - i__ + 1; 22175 sgelq2_(&i__2, &i__1, &a_ref(i__, i__), lda, &tau[i__], &work[1], & 22176 iinfo); 22177 } 22178 22179 work[1] = (real) iws; 22180 return 0; 22181 22182 /* End of SGELQF */ 22183 22184 } /* sgelqf_ */
| int sgemm_ | ( | const char * | transa, | |
| const char * | transb, | |||
| integer * | m, | |||
| integer * | n, | |||
| integer * | k, | |||
| real * | alpha, | |||
| real * | a, | |||
| integer * | lda, | |||
| real * | b, | |||
| integer * | ldb, | |||
| real * | beta, | |||
| real * | c__, | |||
| integer * | ldc | |||
| ) |
Definition at line 1312 of file lapackblas.cpp.
References a_ref, b_ref, c___ref, f2cmax, integer, lsame_(), and xerbla_().
Referenced by sgebrd_(), sgesvd_(), slaed0_(), slaed3_(), slaed7_(), slarfb_(), and sstedc_().
01315 { 01316 /* System generated locals */ 01317 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2, 01318 i__3; 01319 /* Local variables */ 01320 static integer info; 01321 static logical nota, notb; 01322 static real temp; 01323 static integer i__, j, l, ncola; 01324 extern logical lsame_(const char *, const char *); 01325 static integer nrowa, nrowb; 01326 extern /* Subroutine */ int xerbla_(const char *, integer *); 01327 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 01328 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1] 01329 #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1] 01330 /* Purpose 01331 ======= 01332 SGEMM performs one of the matrix-matrix operations 01333 C := alpha*op( A )*op( B ) + beta*C, 01334 where op( X ) is one of 01335 op( X ) = X or op( X ) = X', 01336 alpha and beta are scalars, and A, B and C are matrices, with op( A ) 01337 an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. 01338 Parameters 01339 ========== 01340 TRANSA - CHARACTER*1. 01341 On entry, TRANSA specifies the form of op( A ) to be used in 01342 the matrix multiplication as follows: 01343 TRANSA = 'N' or 'n', op( A ) = A. 01344 TRANSA = 'T' or 't', op( A ) = A'. 01345 TRANSA = 'C' or 'c', op( A ) = A'. 01346 Unchanged on exit. 01347 TRANSB - CHARACTER*1. 01348 On entry, TRANSB specifies the form of op( B ) to be used in 01349 the matrix multiplication as follows: 01350 TRANSB = 'N' or 'n', op( B ) = B. 01351 TRANSB = 'T' or 't', op( B ) = B'. 01352 TRANSB = 'C' or 'c', op( B ) = B'. 01353 Unchanged on exit. 01354 M - INTEGER. 01355 On entry, M specifies the number of rows of the matrix 01356 op( A ) and of the matrix C. M must be at least zero. 01357 Unchanged on exit. 01358 N - INTEGER. 01359 On entry, N specifies the number of columns of the matrix 01360 op( B ) and the number of columns of the matrix C. N must be 01361 at least zero. 01362 Unchanged on exit. 01363 K - INTEGER. 01364 On entry, K specifies the number of columns of the matrix 01365 op( A ) and the number of rows of the matrix op( B ). K must 01366 be at least zero. 01367 Unchanged on exit. 01368 ALPHA - REAL . 01369 On entry, ALPHA specifies the scalar alpha. 01370 Unchanged on exit. 01371 A - REAL array of DIMENSION ( LDA, ka ), where ka is 01372 k when TRANSA = 'N' or 'n', and is m otherwise. 01373 Before entry with TRANSA = 'N' or 'n', the leading m by k 01374 part of the array A must contain the matrix A, otherwise 01375 the leading k by m part of the array A must contain the 01376 matrix A. 01377 Unchanged on exit. 01378 LDA - INTEGER. 01379 On entry, LDA specifies the first dimension of A as declared 01380 in the calling (sub) program. When TRANSA = 'N' or 'n' then 01381 LDA must be at least f2cmax( 1, m ), otherwise LDA must be at 01382 least f2cmax( 1, k ). 01383 Unchanged on exit. 01384 B - REAL array of DIMENSION ( LDB, kb ), where kb is 01385 n when TRANSB = 'N' or 'n', and is k otherwise. 01386 Before entry with TRANSB = 'N' or 'n', the leading k by n 01387 part of the array B must contain the matrix B, otherwise 01388 the leading n by k part of the array B must contain the 01389 matrix B. 01390 Unchanged on exit. 01391 LDB - INTEGER. 01392 On entry, LDB specifies the first dimension of B as declared 01393 in the calling (sub) program. When TRANSB = 'N' or 'n' then 01394 LDB must be at least f2cmax( 1, k ), otherwise LDB must be at 01395 least f2cmax( 1, n ). 01396 Unchanged on exit. 01397 BETA - REAL . 01398 On entry, BETA specifies the scalar beta. When BETA is 01399 supplied as zero then C need not be set on input. 01400 Unchanged on exit. 01401 C - REAL array of DIMENSION ( LDC, n ). 01402 Before entry, the leading m by n part of the array C must 01403 contain the matrix C, except when beta is zero, in which 01404 case C need not be set on entry. 01405 On exit, the array C is overwritten by the m by n matrix 01406 ( alpha*op( A )*op( B ) + beta*C ). 01407 LDC - INTEGER. 01408 On entry, LDC specifies the first dimension of C as declared 01409 in the calling (sub) program. LDC must be at least 01410 f2cmax( 1, m ). 01411 Unchanged on exit. 01412 Level 3 Blas routine. 01413 -- Written on 8-February-1989. 01414 Jack Dongarra, Argonne National Laboratory. 01415 Iain Duff, AERE Harwell. 01416 Jeremy Du Croz, Numerical Algorithms Group Ltd. 01417 Sven Hammarling, Numerical Algorithms Group Ltd. 01418 Set NOTA and NOTB as true if A and B respectively are not 01419 transposed and set NROWA, NCOLA and NROWB as the number of rows 01420 and columns of A and the number of rows of B respectively. 01421 Parameter adjustments */ 01422 a_dim1 = *lda; 01423 a_offset = 1 + a_dim1 * 1; 01424 a -= a_offset; 01425 b_dim1 = *ldb; 01426 b_offset = 1 + b_dim1 * 1; 01427 b -= b_offset; 01428 c_dim1 = *ldc; 01429 c_offset = 1 + c_dim1 * 1; 01430 c__ -= c_offset; 01431 /* Function Body */ 01432 nota = lsame_(transa, "N"); 01433 notb = lsame_(transb, "N"); 01434 if (nota) { 01435 nrowa = *m; 01436 ncola = *k; 01437 } else { 01438 nrowa = *k; 01439 ncola = *m; 01440 } 01441 if (notb) { 01442 nrowb = *k; 01443 } else { 01444 nrowb = *n; 01445 } 01446 /* Test the input parameters. */ 01447 info = 0; 01448 if (! nota && ! lsame_(transa, "C") && ! lsame_( 01449 transa, "T")) { 01450 info = 1; 01451 } else if (! notb && ! lsame_(transb, "C") && ! 01452 lsame_(transb, "T")) { 01453 info = 2; 01454 } else if (*m < 0) { 01455 info = 3; 01456 } else if (*n < 0) { 01457 info = 4; 01458 } else if (*k < 0) { 01459 info = 5; 01460 } else if (*lda < f2cmax(1,nrowa)) { 01461 info = 8; 01462 } else if (*ldb < f2cmax(1,nrowb)) { 01463 info = 10; 01464 } else if (*ldc < f2cmax(1,*m)) { 01465 info = 13; 01466 } 01467 if (info != 0) { 01468 xerbla_("SGEMM ", &info); 01469 return 0; 01470 } 01471 /* Quick return if possible. */ 01472 if (*m == 0 || *n == 0 || (*alpha == 0.f || *k == 0) && *beta == 1.f) { 01473 return 0; 01474 } 01475 /* And if alpha.eq.zero. */ 01476 if (*alpha == 0.f) { 01477 if (*beta == 0.f) { 01478 i__1 = *n; 01479 for (j = 1; j <= i__1; ++j) { 01480 i__2 = *m; 01481 for (i__ = 1; i__ <= i__2; ++i__) { 01482 c___ref(i__, j) = 0.f; 01483 /* L10: */ 01484 } 01485 /* L20: */ 01486 } 01487 } else { 01488 i__1 = *n; 01489 for (j = 1; j <= i__1; ++j) { 01490 i__2 = *m; 01491 for (i__ = 1; i__ <= i__2; ++i__) { 01492 c___ref(i__, j) = *beta * c___ref(i__, j); 01493 /* L30: */ 01494 } 01495 /* L40: */ 01496 } 01497 } 01498 return 0; 01499 } 01500 /* Start the operations. */ 01501 if (notb) { 01502 if (nota) { 01503 /* Form C := alpha*A*B + beta*C. */ 01504 i__1 = *n; 01505 for (j = 1; j <= i__1; ++j) { 01506 if (*beta == 0.f) { 01507 i__2 = *m; 01508 for (i__ = 1; i__ <= i__2; ++i__) { 01509 c___ref(i__, j) = 0.f; 01510 /* L50: */ 01511 } 01512 } else if (*beta != 1.f) { 01513 i__2 = *m; 01514 for (i__ = 1; i__ <= i__2; ++i__) { 01515 c___ref(i__, j) = *beta * c___ref(i__, j); 01516 /* L60: */ 01517 } 01518 } 01519 i__2 = *k; 01520 for (l = 1; l <= i__2; ++l) { 01521 if (b_ref(l, j) != 0.f) { 01522 temp = *alpha * b_ref(l, j); 01523 i__3 = *m; 01524 for (i__ = 1; i__ <= i__3; ++i__) { 01525 c___ref(i__, j) = c___ref(i__, j) + temp * a_ref( 01526 i__, l); 01527 /* L70: */ 01528 } 01529 } 01530 /* L80: */ 01531 } 01532 /* L90: */ 01533 } 01534 } else { 01535 /* Form C := alpha*A'*B + beta*C */ 01536 i__1 = *n; 01537 for (j = 1; j <= i__1; ++j) { 01538 i__2 = *m; 01539 for (i__ = 1; i__ <= i__2; ++i__) { 01540 temp = 0.f; 01541 i__3 = *k; 01542 for (l = 1; l <= i__3; ++l) { 01543 temp += a_ref(l, i__) * b_ref(l, j); 01544 /* L100: */ 01545 } 01546 if (*beta == 0.f) { 01547 c___ref(i__, j) = *alpha * temp; 01548 } else { 01549 c___ref(i__, j) = *alpha * temp + *beta * c___ref(i__, 01550 j); 01551 } 01552 /* L110: */ 01553 } 01554 /* L120: */ 01555 } 01556 } 01557 } else { 01558 if (nota) { 01559 /* Form C := alpha*A*B' + beta*C */ 01560 i__1 = *n; 01561 for (j = 1; j <= i__1; ++j) { 01562 if (*beta == 0.f) { 01563 i__2 = *m; 01564 for (i__ = 1; i__ <= i__2; ++i__) { 01565 c___ref(i__, j) = 0.f; 01566 /* L130: */ 01567 } 01568 } else if (*beta != 1.f) { 01569 i__2 = *m; 01570 for (i__ = 1; i__ <= i__2; ++i__) { 01571 c___ref(i__, j) = *beta * c___ref(i__, j); 01572 /* L140: */ 01573 } 01574 } 01575 i__2 = *k; 01576 for (l = 1; l <= i__2; ++l) { 01577 if (b_ref(j, l) != 0.f) { 01578 temp = *alpha * b_ref(j, l); 01579 i__3 = *m; 01580 for (i__ = 1; i__ <= i__3; ++i__) { 01581 c___ref(i__, j) = c___ref(i__, j) + temp * a_ref( 01582 i__, l); 01583 /* L150: */ 01584 } 01585 } 01586 /* L160: */ 01587 } 01588 /* L170: */ 01589 } 01590 } else { 01591 /* Form C := alpha*A'*B' + beta*C */ 01592 i__1 = *n; 01593 for (j = 1; j <= i__1; ++j) { 01594 i__2 = *m; 01595 for (i__ = 1; i__ <= i__2; ++i__) { 01596 temp = 0.f; 01597 i__3 = *k; 01598 for (l = 1; l <= i__3; ++l) { 01599 temp += a_ref(l, i__) * b_ref(j, l); 01600 /* L180: */ 01601 } 01602 if (*beta == 0.f) { 01603 c___ref(i__, j) = *alpha * temp; 01604 } else { 01605 c___ref(i__, j) = *alpha * temp + *beta * c___ref(i__, 01606 j); 01607 } 01608 /* L190: */ 01609 } 01610 /* L200: */ 01611 } 01612 } 01613 } 01614 return 0; 01615 /* End of SGEMM . */ 01616 } /* sgemm_ */
| int sgemv_ | ( | const char * | trans, | |
| integer * | m, | |||
| integer * | n, | |||
| real * | alpha, | |||
| real * | a, | |||
| integer * | lda, | |||
| real * | x, | |||
| integer * | incx, | |||
| real * | beta, | |||
| real * | y, | |||
| integer * | incy | |||
| ) |
Definition at line 1624 of file lapackblas.cpp.
References a_ref, f2cmax, integer, lsame_(), and xerbla_().
Referenced by EMAN::PCAlarge::analyze(), EMAN::PCA::dopca_lan(), EMAN::PCA::Lanczos(), EMAN::PCAlarge::Lanczos(), slabrd_(), slaeda_(), slarf_(), slarft_(), and slatrd_().
01627 { 01628 /* System generated locals */ 01629 integer a_dim1, a_offset, i__1, i__2; 01630 /* Local variables */ 01631 static integer info; 01632 static real temp; 01633 static integer lenx, leny, i__, j; 01634 extern logical lsame_(const char *, const char *); 01635 static integer ix, iy, jx, jy, kx, ky; 01636 extern /* Subroutine */ int xerbla_(const char *, integer *); 01637 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 01638 /* Purpose 01639 ======= 01640 SGEMV performs one of the matrix-vector operations 01641 y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, 01642 where alpha and beta are scalars, x and y are vectors and A is an 01643 m by n matrix. 01644 Parameters 01645 ========== 01646 TRANS - CHARACTER*1. 01647 On entry, TRANS specifies the operation to be performed as 01648 follows: 01649 TRANS = 'N' or 'n' y := alpha*A*x + beta*y. 01650 TRANS = 'T' or 't' y := alpha*A'*x + beta*y. 01651 TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. 01652 Unchanged on exit. 01653 M - INTEGER. 01654 On entry, M specifies the number of rows of the matrix A. 01655 M must be at least zero. 01656 Unchanged on exit. 01657 N - INTEGER. 01658 On entry, N specifies the number of columns of the matrix A. 01659 N must be at least zero. 01660 Unchanged on exit. 01661 ALPHA - REAL . 01662 On entry, ALPHA specifies the scalar alpha. 01663 Unchanged on exit. 01664 A - REAL array of DIMENSION ( LDA, n ). 01665 Before entry, the leading m by n part of the array A must 01666 contain the matrix of coefficients. 01667 Unchanged on exit. 01668 LDA - INTEGER. 01669 On entry, LDA specifies the first dimension of A as declared 01670 in the calling (sub) program. LDA must be at least 01671 f2cmax( 1, m ). 01672 Unchanged on exit. 01673 X - REAL array of DIMENSION at least 01674 ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' 01675 and at least 01676 ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. 01677 Before entry, the incremented array X must contain the 01678 vector x. 01679 Unchanged on exit. 01680 INCX - INTEGER. 01681 On entry, INCX specifies the increment for the elements of 01682 X. INCX must not be zero. 01683 Unchanged on exit. 01684 BETA - REAL . 01685 On entry, BETA specifies the scalar beta. When BETA is 01686 supplied as zero then Y need not be set on input. 01687 Unchanged on exit. 01688 Y - REAL array of DIMENSION at least 01689 ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' 01690 and at least 01691 ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. 01692 Before entry with BETA non-zero, the incremented array Y 01693 must contain the vector y. On exit, Y is overwritten by the 01694 updated vector y. 01695 INCY - INTEGER. 01696 On entry, INCY specifies the increment for the elements of 01697 Y. INCY must not be zero. 01698 Unchanged on exit. 01699 Level 2 Blas routine. 01700 -- Written on 22-October-1986. 01701 Jack Dongarra, Argonne National Lab. 01702 Jeremy Du Croz, Nag Central Office. 01703 Sven Hammarling, Nag Central Office. 01704 Richard Hanson, Sandia National Labs. 01705 Test the input parameters. 01706 Parameter adjustments */ 01707 a_dim1 = *lda; 01708 a_offset = 1 + a_dim1 * 1; 01709 a -= a_offset; 01710 --x; 01711 --y; 01712 /* Function Body */ 01713 info = 0; 01714 if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C") 01715 ) { 01716 info = 1; 01717 } else if (*m < 0) { 01718 info = 2; 01719 } else if (*n < 0) { 01720 info = 3; 01721 } else if (*lda < f2cmax(1,*m)) { 01722 info = 6; 01723 } else if (*incx == 0) { 01724 info = 8; 01725 } else if (*incy == 0) { 01726 info = 11; 01727 } 01728 if (info != 0) { 01729 xerbla_("SGEMV ", &info); 01730 return 0; 01731 } 01732 /* Quick return if possible. */ 01733 if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f) { 01734 return 0; 01735 } 01736 /* Set LENX and LENY, the lengths of the vectors x and y, and set 01737 up the start points in X and Y. */ 01738 if (lsame_(trans, "N")) { 01739 lenx = *n; 01740 leny = *m; 01741 } else { 01742 lenx = *m; 01743 leny = *n; 01744 } 01745 if (*incx > 0) { 01746 kx = 1; 01747 } else { 01748 kx = 1 - (lenx - 1) * *incx; 01749 } 01750 if (*incy > 0) { 01751 ky = 1; 01752 } else { 01753 ky = 1 - (leny - 1) * *incy; 01754 } 01755 /* Start the operations. In this version the elements of A are 01756 accessed sequentially with one pass through A. 01757 First form y := beta*y. */ 01758 if (*beta != 1.f) { 01759 if (*incy == 1) { 01760 if (*beta == 0.f) { 01761 i__1 = leny; 01762 for (i__ = 1; i__ <= i__1; ++i__) { 01763 y[i__] = 0.f; 01764 /* L10: */ 01765 } 01766 } else { 01767 i__1 = leny; 01768 for (i__ = 1; i__ <= i__1; ++i__) { 01769 y[i__] = *beta * y[i__]; 01770 /* L20: */ 01771 } 01772 } 01773 } else { 01774 iy = ky; 01775 if (*beta == 0.f) { 01776 i__1 = leny; 01777 for (i__ = 1; i__ <= i__1; ++i__) { 01778 y[iy] = 0.f; 01779 iy += *incy; 01780 /* L30: */ 01781 } 01782 } else { 01783 i__1 = leny; 01784 for (i__ = 1; i__ <= i__1; ++i__) { 01785 y[iy] = *beta * y[iy]; 01786 iy += *incy; 01787 /* L40: */ 01788 } 01789 } 01790 } 01791 } 01792 if (*alpha == 0.f) { 01793 return 0; 01794 } 01795 if (lsame_(trans, "N")) { 01796 /* Form y := alpha*A*x + y. */ 01797 jx = kx; 01798 if (*incy == 1) { 01799 i__1 = *n; 01800 for (j = 1; j <= i__1; ++j) { 01801 if (x[jx] != 0.f) { 01802 temp = *alpha * x[jx]; 01803 i__2 = *m; 01804 for (i__ = 1; i__ <= i__2; ++i__) { 01805 y[i__] += temp * a_ref(i__, j); 01806 /* L50: */ 01807 } 01808 } 01809 jx += *incx; 01810 /* L60: */ 01811 } 01812 } else { 01813 i__1 = *n; 01814 for (j = 1; j <= i__1; ++j) { 01815 if (x[jx] != 0.f) { 01816 temp = *alpha * x[jx]; 01817 iy = ky; 01818 i__2 = *m; 01819 for (i__ = 1; i__ <= i__2; ++i__) { 01820 y[iy] += temp * a_ref(i__, j); 01821 iy += *incy; 01822 /* L70: */ 01823 } 01824 } 01825 jx += *incx; 01826 /* L80: */ 01827 } 01828 } 01829 } else { 01830 /* Form y := alpha*A'*x + y. */ 01831 jy = ky; 01832 if (*incx == 1) { 01833 i__1 = *n; 01834 for (j = 1; j <= i__1; ++j) { 01835 temp = 0.f; 01836 i__2 = *m; 01837 for (i__ = 1; i__ <= i__2; ++i__) { 01838 temp += a_ref(i__, j) * x[i__]; 01839 /* L90: */ 01840 } 01841 y[jy] += *alpha * temp; 01842 jy += *incy; 01843 /* L100: */ 01844 } 01845 } else { 01846 i__1 = *n; 01847 for (j = 1; j <= i__1; ++j) { 01848 temp = 0.f; 01849 ix = kx; 01850 i__2 = *m; 01851 for (i__ = 1; i__ <= i__2; ++i__) { 01852 temp += a_ref(i__, j) * x[ix]; 01853 ix += *incx; 01854 /* L110: */ 01855 } 01856 y[jy] += *alpha * temp; 01857 jy += *incy; 01858 /* L120: */ 01859 } 01860 } 01861 } 01862 return 0; 01863 /* End of SGEMV . */ 01864 } /* sgemv_ */
| int sgeqr2_ | ( | integer * | m, | |
| integer * | n, | |||
| real * | a, | |||
| integer * | lda, | |||
| real * | tau, | |||
| real * | work, | |||
| integer * | info | |||
| ) |
Definition at line 24684 of file lapackblas.cpp.
References a_ref, f2cmax, f2cmin, integer, slarf_(), slarfg_(), and xerbla_().
Referenced by sgeqrf_().
24686 { 24687 /* -- LAPACK routine (version 3.0) -- 24688 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 24689 Courant Institute, Argonne National Lab, and Rice University 24690 February 29, 1992 24691 24692 24693 Purpose 24694 ======= 24695 24696 SGEQR2 computes a QR factorization of a real m by n matrix A: 24697 A = Q * R. 24698 24699 Arguments 24700 ========= 24701 24702 M (input) INTEGER 24703 The number of rows of the matrix A. M >= 0. 24704 24705 N (input) INTEGER 24706 The number of columns of the matrix A. N >= 0. 24707 24708 A (input/output) REAL array, dimension (LDA,N) 24709 On entry, the m by n matrix A. 24710 On exit, the elements on and above the diagonal of the array 24711 contain the min(m,n) by n upper trapezoidal matrix R (R is 24712 upper triangular if m >= n); the elements below the diagonal, 24713 with the array TAU, represent the orthogonal matrix Q as a 24714 product of elementary reflectors (see Further Details). 24715