Deprecate ?GELQS and ?GEQRS from TESTING/LIN (Reference-LAPACK PR 900) (#4307)
* Move ?GELQS and ?GEQRS from TESTING/LIN to DEPRECATED (Reference-LAPACK PR 900) * Add f2c-converted versions of ?GELQS and ?GEQRS
This commit is contained in:
parent
00ef1bb58a
commit
58427ff74d
|
@ -438,15 +438,19 @@ endif()
|
|||
|
||||
if(BUILD_LAPACK_DEPRECATED)
|
||||
list(APPEND SLASRC DEPRECATED/sgegs.f DEPRECATED/sgegv.f
|
||||
DEPRECATED/sgelqs.f DEPRECATED/sgeqrs.f
|
||||
DEPRECATED/sgeqpf.f DEPRECATED/sgelsx.f DEPRECATED/sggsvd.f
|
||||
DEPRECATED/sggsvp.f DEPRECATED/slahrd.f DEPRECATED/slatzm.f DEPRECATED/stzrqf.f)
|
||||
list(APPEND DLASRC DEPRECATED/dgegs.f DEPRECATED/dgegv.f
|
||||
DEPRECATED/dgelqs.f DEPRECATED/dgeqrs.f
|
||||
DEPRECATED/dgeqpf.f DEPRECATED/dgelsx.f DEPRECATED/dggsvd.f
|
||||
DEPRECATED/dggsvp.f DEPRECATED/dlahrd.f DEPRECATED/dlatzm.f DEPRECATED/dtzrqf.f)
|
||||
list(APPEND CLASRC DEPRECATED/cgegs.f DEPRECATED/cgegv.f
|
||||
DEPRECATED/cgelqs.f DEPRECATED/cgeqrs.f
|
||||
DEPRECATED/cgeqpf.f DEPRECATED/cgelsx.f DEPRECATED/cggsvd.f
|
||||
DEPRECATED/cggsvp.f DEPRECATED/clahrd.f DEPRECATED/clatzm.f DEPRECATED/ctzrqf.f)
|
||||
list(APPEND ZLASRC DEPRECATED/zgegs.f DEPRECATED/zgegv.f
|
||||
DEPRECATED/zgelqs.f DEPRECATED/zgeqrs.f
|
||||
DEPRECATED/zgeqpf.f DEPRECATED/zgelsx.f DEPRECATED/zggsvd.f
|
||||
DEPRECATED/zggsvp.f DEPRECATED/zlahrd.f DEPRECATED/zlatzm.f DEPRECATED/ztzrqf.f)
|
||||
message(STATUS "Building deprecated routines")
|
||||
|
@ -935,15 +939,19 @@ endif()
|
|||
|
||||
if(BUILD_LAPACK_DEPRECATED)
|
||||
list(APPEND SLASRC DEPRECATED/sgegs.c DEPRECATED/sgegv.c
|
||||
DEPRECATED/sgelqs.c DEPRECATED/sgeqrs.c
|
||||
DEPRECATED/sgeqpf.c DEPRECATED/sgelsx.c DEPRECATED/sggsvd.c
|
||||
DEPRECATED/sggsvp.c DEPRECATED/slahrd.c DEPRECATED/slatzm.c DEPRECATED/stzrqf.c)
|
||||
list(APPEND DLASRC DEPRECATED/dgegs.c DEPRECATED/dgegv.c
|
||||
DEPRECATED/dgelqs.c DEPRECATED/dgeqrs.c
|
||||
DEPRECATED/dgeqpf.c DEPRECATED/dgelsx.c DEPRECATED/dggsvd.c
|
||||
DEPRECATED/dggsvp.c DEPRECATED/dlahrd.c DEPRECATED/dlatzm.c DEPRECATED/dtzrqf.c)
|
||||
list(APPEND CLASRC DEPRECATED/cgegs.c DEPRECATED/cgegv.c
|
||||
DEPRECATED/cgelqs.c DEPRECATED/cgeqrs.c
|
||||
DEPRECATED/cgeqpf.c DEPRECATED/cgelsx.c DEPRECATED/cggsvd.c
|
||||
DEPRECATED/cggsvp.c DEPRECATED/clahrd.c DEPRECATED/clatzm.c DEPRECATED/ctzrqf.c)
|
||||
list(APPEND ZLASRC DEPRECATED/zgegs.c DEPRECATED/zgegv.c
|
||||
DEPRECATED/zgelqs.c DEPRECATED/zgeqrs.c
|
||||
DEPRECATED/zgeqpf.c DEPRECATED/zgelsx.c DEPRECATED/zggsvd.c
|
||||
DEPRECATED/zggsvp.c DEPRECATED/zlahrd.c DEPRECATED/zlatzm.c DEPRECATED/ztzrqf.c)
|
||||
message(STATUS "Building deprecated routines")
|
||||
|
|
|
@ -0,0 +1,479 @@
|
|||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <stdio.h>
|
||||
#include <complex.h>
|
||||
#ifdef complex
|
||||
#undef complex
|
||||
#endif
|
||||
#ifdef I
|
||||
#undef I
|
||||
#endif
|
||||
|
||||
#if defined(_WIN64)
|
||||
typedef long long BLASLONG;
|
||||
typedef unsigned long long BLASULONG;
|
||||
#else
|
||||
typedef long BLASLONG;
|
||||
typedef unsigned long BLASULONG;
|
||||
#endif
|
||||
|
||||
#ifdef LAPACK_ILP64
|
||||
typedef BLASLONG blasint;
|
||||
#if defined(_WIN64)
|
||||
#define blasabs(x) llabs(x)
|
||||
#else
|
||||
#define blasabs(x) labs(x)
|
||||
#endif
|
||||
#else
|
||||
typedef int blasint;
|
||||
#define blasabs(x) abs(x)
|
||||
#endif
|
||||
|
||||
typedef blasint integer;
|
||||
|
||||
typedef unsigned int uinteger;
|
||||
typedef char *address;
|
||||
typedef short int shortint;
|
||||
typedef float real;
|
||||
typedef double doublereal;
|
||||
typedef struct { real r, i; } complex;
|
||||
typedef struct { doublereal r, i; } doublecomplex;
|
||||
#ifdef _MSC_VER
|
||||
static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
|
||||
static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
|
||||
static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
|
||||
static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
|
||||
#else
|
||||
static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
|
||||
static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
|
||||
#endif
|
||||
#define pCf(z) (*_pCf(z))
|
||||
#define pCd(z) (*_pCd(z))
|
||||
typedef int logical;
|
||||
typedef short int shortlogical;
|
||||
typedef char logical1;
|
||||
typedef char integer1;
|
||||
|
||||
#define TRUE_ (1)
|
||||
#define FALSE_ (0)
|
||||
|
||||
/* Extern is for use with -E */
|
||||
#ifndef Extern
|
||||
#define Extern extern
|
||||
#endif
|
||||
|
||||
/* I/O stuff */
|
||||
|
||||
typedef int flag;
|
||||
typedef int ftnlen;
|
||||
typedef int ftnint;
|
||||
|
||||
/*external read, write*/
|
||||
typedef struct
|
||||
{ flag cierr;
|
||||
ftnint ciunit;
|
||||
flag ciend;
|
||||
char *cifmt;
|
||||
ftnint cirec;
|
||||
} cilist;
|
||||
|
||||
/*internal read, write*/
|
||||
typedef struct
|
||||
{ flag icierr;
|
||||
char *iciunit;
|
||||
flag iciend;
|
||||
char *icifmt;
|
||||
ftnint icirlen;
|
||||
ftnint icirnum;
|
||||
} icilist;
|
||||
|
||||
/*open*/
|
||||
typedef struct
|
||||
{ flag oerr;
|
||||
ftnint ounit;
|
||||
char *ofnm;
|
||||
ftnlen ofnmlen;
|
||||
char *osta;
|
||||
char *oacc;
|
||||
char *ofm;
|
||||
ftnint orl;
|
||||
char *oblnk;
|
||||
} olist;
|
||||
|
||||
/*close*/
|
||||
typedef struct
|
||||
{ flag cerr;
|
||||
ftnint cunit;
|
||||
char *csta;
|
||||
} cllist;
|
||||
|
||||
/*rewind, backspace, endfile*/
|
||||
typedef struct
|
||||
{ flag aerr;
|
||||
ftnint aunit;
|
||||
} alist;
|
||||
|
||||
/* inquire */
|
||||
typedef struct
|
||||
{ flag inerr;
|
||||
ftnint inunit;
|
||||
char *infile;
|
||||
ftnlen infilen;
|
||||
ftnint *inex; /*parameters in standard's order*/
|
||||
ftnint *inopen;
|
||||
ftnint *innum;
|
||||
ftnint *innamed;
|
||||
char *inname;
|
||||
ftnlen innamlen;
|
||||
char *inacc;
|
||||
ftnlen inacclen;
|
||||
char *inseq;
|
||||
ftnlen inseqlen;
|
||||
char *indir;
|
||||
ftnlen indirlen;
|
||||
char *infmt;
|
||||
ftnlen infmtlen;
|
||||
char *inform;
|
||||
ftnint informlen;
|
||||
char *inunf;
|
||||
ftnlen inunflen;
|
||||
ftnint *inrecl;
|
||||
ftnint *innrec;
|
||||
char *inblank;
|
||||
ftnlen inblanklen;
|
||||
} inlist;
|
||||
|
||||
#define VOID void
|
||||
|
||||
union Multitype { /* for multiple entry points */
|
||||
integer1 g;
|
||||
shortint h;
|
||||
integer i;
|
||||
/* longint j; */
|
||||
real r;
|
||||
doublereal d;
|
||||
complex c;
|
||||
doublecomplex z;
|
||||
};
|
||||
|
||||
typedef union Multitype Multitype;
|
||||
|
||||
struct Vardesc { /* for Namelist */
|
||||
char *name;
|
||||
char *addr;
|
||||
ftnlen *dims;
|
||||
int type;
|
||||
};
|
||||
typedef struct Vardesc Vardesc;
|
||||
|
||||
struct Namelist {
|
||||
char *name;
|
||||
Vardesc **vars;
|
||||
int nvars;
|
||||
};
|
||||
typedef struct Namelist Namelist;
|
||||
|
||||
#define abs(x) ((x) >= 0 ? (x) : -(x))
|
||||
#define dabs(x) (fabs(x))
|
||||
#define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
|
||||
#define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
|
||||
#define dmin(a,b) (f2cmin(a,b))
|
||||
#define dmax(a,b) (f2cmax(a,b))
|
||||
#define bit_test(a,b) ((a) >> (b) & 1)
|
||||
#define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
|
||||
#define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
|
||||
|
||||
#define abort_() { sig_die("Fortran abort routine called", 1); }
|
||||
#define c_abs(z) (cabsf(Cf(z)))
|
||||
#define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
|
||||
#ifdef _MSC_VER
|
||||
#define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
|
||||
#define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/Cd(b)._Val[1]);}
|
||||
#else
|
||||
#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
|
||||
#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
|
||||
#endif
|
||||
#define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
|
||||
#define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
|
||||
#define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
|
||||
//#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
|
||||
#define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
|
||||
#define d_abs(x) (fabs(*(x)))
|
||||
#define d_acos(x) (acos(*(x)))
|
||||
#define d_asin(x) (asin(*(x)))
|
||||
#define d_atan(x) (atan(*(x)))
|
||||
#define d_atn2(x, y) (atan2(*(x),*(y)))
|
||||
#define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
|
||||
#define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
|
||||
#define d_cos(x) (cos(*(x)))
|
||||
#define d_cosh(x) (cosh(*(x)))
|
||||
#define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
|
||||
#define d_exp(x) (exp(*(x)))
|
||||
#define d_imag(z) (cimag(Cd(z)))
|
||||
#define r_imag(z) (cimagf(Cf(z)))
|
||||
#define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define d_log(x) (log(*(x)))
|
||||
#define d_mod(x, y) (fmod(*(x), *(y)))
|
||||
#define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
|
||||
#define d_nint(x) u_nint(*(x))
|
||||
#define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
|
||||
#define d_sign(a,b) u_sign(*(a),*(b))
|
||||
#define r_sign(a,b) u_sign(*(a),*(b))
|
||||
#define d_sin(x) (sin(*(x)))
|
||||
#define d_sinh(x) (sinh(*(x)))
|
||||
#define d_sqrt(x) (sqrt(*(x)))
|
||||
#define d_tan(x) (tan(*(x)))
|
||||
#define d_tanh(x) (tanh(*(x)))
|
||||
#define i_abs(x) abs(*(x))
|
||||
#define i_dnnt(x) ((integer)u_nint(*(x)))
|
||||
#define i_len(s, n) (n)
|
||||
#define i_nint(x) ((integer)u_nint(*(x)))
|
||||
#define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
|
||||
#define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
|
||||
#define pow_si(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_ri(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_di(B,E) dpow_ui(*(B),*(E))
|
||||
#define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
|
||||
#define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
|
||||
#define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
|
||||
#define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
|
||||
#define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
|
||||
#define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
|
||||
#define sig_die(s, kill) { exit(1); }
|
||||
#define s_stop(s, n) {exit(0);}
|
||||
#define z_abs(z) (cabs(Cd(z)))
|
||||
#define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
|
||||
#define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
|
||||
#define myexit_() break;
|
||||
#define mycycle_() continue;
|
||||
#define myceiling_(w) {ceil(w)}
|
||||
#define myhuge_(w) {HUGE_VAL}
|
||||
#define mymaxloc_(w,s,e,n) dmaxloc_(w,*(s),*(e),n)
|
||||
|
||||
/* procedure parameter types for -A and -C++ */
|
||||
|
||||
#define F2C_proc_par_types 1
|
||||
#ifdef __cplusplus
|
||||
typedef logical (*L_fp)(...);
|
||||
#else
|
||||
typedef logical (*L_fp)();
|
||||
#endif
|
||||
|
||||
/* -- translated by f2c (version 20000121).
|
||||
You must link the resulting object file with the libraries:
|
||||
-lf2c -lm (in that order)
|
||||
*/
|
||||
|
||||
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static complex c_b1 = {0.f,0.f};
|
||||
static complex c_b2 = {1.f,0.f};
|
||||
|
||||
/* > \brief \b CGELQS */
|
||||
|
||||
/* =========== DOCUMENTATION =========== */
|
||||
|
||||
/* Online html documentation available at */
|
||||
/* http://www.netlib.org/lapack/explore-html/ */
|
||||
|
||||
/* Definition: */
|
||||
/* =========== */
|
||||
|
||||
/* SUBROUTINE CGELQS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, */
|
||||
/* INFO ) */
|
||||
|
||||
/* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS */
|
||||
/* COMPLEX A( LDA, * ), B( LDB, * ), TAU( * ), */
|
||||
/* $ WORK( LWORK ) */
|
||||
|
||||
|
||||
/* > \par Purpose: */
|
||||
/* ============= */
|
||||
/* > */
|
||||
/* > \verbatim */
|
||||
/* > */
|
||||
/* > Compute a minimum-norm solution */
|
||||
/* > f2cmin || A*X - B || */
|
||||
/* > using the LQ factorization */
|
||||
/* > A = L*Q */
|
||||
/* > computed by CGELQF. */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Arguments: */
|
||||
/* ========== */
|
||||
|
||||
/* > \param[in] M */
|
||||
/* > \verbatim */
|
||||
/* > M is INTEGER */
|
||||
/* > The number of rows of the matrix A. M >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] N */
|
||||
/* > \verbatim */
|
||||
/* > N is INTEGER */
|
||||
/* > The number of columns of the matrix A. N >= M >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] NRHS */
|
||||
/* > \verbatim */
|
||||
/* > NRHS is INTEGER */
|
||||
/* > The number of columns of B. NRHS >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] A */
|
||||
/* > \verbatim */
|
||||
/* > A is COMPLEX array, dimension (LDA,N) */
|
||||
/* > Details of the LQ factorization of the original matrix A as */
|
||||
/* > returned by CGELQF. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDA */
|
||||
/* > \verbatim */
|
||||
/* > LDA is INTEGER */
|
||||
/* > The leading dimension of the array A. LDA >= M. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] TAU */
|
||||
/* > \verbatim */
|
||||
/* > TAU is COMPLEX array, dimension (M) */
|
||||
/* > Details of the orthogonal matrix Q. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in,out] B */
|
||||
/* > \verbatim */
|
||||
/* > B is COMPLEX array, dimension (LDB,NRHS) */
|
||||
/* > On entry, the m-by-nrhs right hand side matrix B. */
|
||||
/* > On exit, the n-by-nrhs solution matrix X. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDB */
|
||||
/* > \verbatim */
|
||||
/* > LDB is INTEGER */
|
||||
/* > The leading dimension of the array B. LDB >= N. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] WORK */
|
||||
/* > \verbatim */
|
||||
/* > WORK is COMPLEX array, dimension (LWORK) */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LWORK */
|
||||
/* > \verbatim */
|
||||
/* > LWORK is INTEGER */
|
||||
/* > The length of the array WORK. LWORK must be at least NRHS, */
|
||||
/* > and should be at least NRHS*NB, where NB is the block size */
|
||||
/* > for this environment. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] INFO */
|
||||
/* > \verbatim */
|
||||
/* > INFO is INTEGER */
|
||||
/* > = 0: successful exit */
|
||||
/* > < 0: if INFO = -i, the i-th argument had an illegal value */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Authors: */
|
||||
/* ======== */
|
||||
|
||||
/* > \author Univ. of Tennessee */
|
||||
/* > \author Univ. of California Berkeley */
|
||||
/* > \author Univ. of Colorado Denver */
|
||||
/* > \author NAG Ltd. */
|
||||
|
||||
/* > \ingroup complex_lin */
|
||||
|
||||
/* ===================================================================== */
|
||||
/* Subroutine */ int cgelqs_(integer *m, integer *n, integer *nrhs, complex *
|
||||
a, integer *lda, complex *tau, complex *b, integer *ldb, complex *
|
||||
work, integer *lwork, integer *info)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, b_dim1, b_offset, i__1;
|
||||
|
||||
/* Local variables */
|
||||
extern /* Subroutine */ int ctrsm_(char *, char *, char *, char *,
|
||||
integer *, integer *, complex *, complex *, integer *, complex *,
|
||||
integer *), claset_(char *,
|
||||
integer *, integer *, complex *, complex *, complex *, integer *), xerbla_(char *, integer *), cunmlq_(char *, char
|
||||
*, integer *, integer *, integer *, complex *, integer *, complex
|
||||
*, complex *, integer *, complex *, integer *, integer *);
|
||||
|
||||
|
||||
/* -- LAPACK test routine -- */
|
||||
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
|
||||
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
|
||||
|
||||
|
||||
/* ===================================================================== */
|
||||
|
||||
|
||||
/* Test the input parameters. */
|
||||
|
||||
/* Parameter adjustments */
|
||||
a_dim1 = *lda;
|
||||
a_offset = 1 + a_dim1 * 1;
|
||||
a -= a_offset;
|
||||
--tau;
|
||||
b_dim1 = *ldb;
|
||||
b_offset = 1 + b_dim1 * 1;
|
||||
b -= b_offset;
|
||||
--work;
|
||||
|
||||
/* Function Body */
|
||||
*info = 0;
|
||||
if (*m < 0) {
|
||||
*info = -1;
|
||||
} else if (*n < 0 || *m > *n) {
|
||||
*info = -2;
|
||||
} else if (*nrhs < 0) {
|
||||
*info = -3;
|
||||
} else if (*lda < f2cmax(1,*m)) {
|
||||
*info = -5;
|
||||
} else if (*ldb < f2cmax(1,*n)) {
|
||||
*info = -8;
|
||||
} else if (*lwork < 1 || *lwork < *nrhs && *m > 0 && *n > 0) {
|
||||
*info = -10;
|
||||
}
|
||||
if (*info != 0) {
|
||||
i__1 = -(*info);
|
||||
xerbla_("CGELQS", &i__1);
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Quick return if possible */
|
||||
|
||||
if (*n == 0 || *nrhs == 0 || *m == 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Solve L*X = B(1:m,:) */
|
||||
|
||||
ctrsm_("Left", "Lower", "No transpose", "Non-unit", m, nrhs, &c_b2, &a[
|
||||
a_offset], lda, &b[b_offset], ldb);
|
||||
|
||||
/* Set B(m+1:n,:) to zero */
|
||||
|
||||
if (*m < *n) {
|
||||
i__1 = *n - *m;
|
||||
claset_("Full", &i__1, nrhs, &c_b1, &c_b1, &b[*m + 1 + b_dim1], ldb);
|
||||
}
|
||||
|
||||
/* B := Q' * B */
|
||||
|
||||
cunmlq_("Left", "Conjugate transpose", n, nrhs, m, &a[a_offset], lda, &
|
||||
tau[1], &b[b_offset], ldb, &work[1], lwork, info);
|
||||
|
||||
return 0;
|
||||
|
||||
/* End of CGELQS */
|
||||
|
||||
} /* cgelqs_ */
|
||||
|
|
@ -0,0 +1,471 @@
|
|||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <stdio.h>
|
||||
#include <complex.h>
|
||||
#ifdef complex
|
||||
#undef complex
|
||||
#endif
|
||||
#ifdef I
|
||||
#undef I
|
||||
#endif
|
||||
|
||||
#if defined(_WIN64)
|
||||
typedef long long BLASLONG;
|
||||
typedef unsigned long long BLASULONG;
|
||||
#else
|
||||
typedef long BLASLONG;
|
||||
typedef unsigned long BLASULONG;
|
||||
#endif
|
||||
|
||||
#ifdef LAPACK_ILP64
|
||||
typedef BLASLONG blasint;
|
||||
#if defined(_WIN64)
|
||||
#define blasabs(x) llabs(x)
|
||||
#else
|
||||
#define blasabs(x) labs(x)
|
||||
#endif
|
||||
#else
|
||||
typedef int blasint;
|
||||
#define blasabs(x) abs(x)
|
||||
#endif
|
||||
|
||||
typedef blasint integer;
|
||||
|
||||
typedef unsigned int uinteger;
|
||||
typedef char *address;
|
||||
typedef short int shortint;
|
||||
typedef float real;
|
||||
typedef double doublereal;
|
||||
typedef struct { real r, i; } complex;
|
||||
typedef struct { doublereal r, i; } doublecomplex;
|
||||
#ifdef _MSC_VER
|
||||
static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
|
||||
static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
|
||||
static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
|
||||
static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
|
||||
#else
|
||||
static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
|
||||
static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
|
||||
#endif
|
||||
#define pCf(z) (*_pCf(z))
|
||||
#define pCd(z) (*_pCd(z))
|
||||
typedef int logical;
|
||||
typedef short int shortlogical;
|
||||
typedef char logical1;
|
||||
typedef char integer1;
|
||||
|
||||
#define TRUE_ (1)
|
||||
#define FALSE_ (0)
|
||||
|
||||
/* Extern is for use with -E */
|
||||
#ifndef Extern
|
||||
#define Extern extern
|
||||
#endif
|
||||
|
||||
/* I/O stuff */
|
||||
|
||||
typedef int flag;
|
||||
typedef int ftnlen;
|
||||
typedef int ftnint;
|
||||
|
||||
/*external read, write*/
|
||||
typedef struct
|
||||
{ flag cierr;
|
||||
ftnint ciunit;
|
||||
flag ciend;
|
||||
char *cifmt;
|
||||
ftnint cirec;
|
||||
} cilist;
|
||||
|
||||
/*internal read, write*/
|
||||
typedef struct
|
||||
{ flag icierr;
|
||||
char *iciunit;
|
||||
flag iciend;
|
||||
char *icifmt;
|
||||
ftnint icirlen;
|
||||
ftnint icirnum;
|
||||
} icilist;
|
||||
|
||||
/*open*/
|
||||
typedef struct
|
||||
{ flag oerr;
|
||||
ftnint ounit;
|
||||
char *ofnm;
|
||||
ftnlen ofnmlen;
|
||||
char *osta;
|
||||
char *oacc;
|
||||
char *ofm;
|
||||
ftnint orl;
|
||||
char *oblnk;
|
||||
} olist;
|
||||
|
||||
/*close*/
|
||||
typedef struct
|
||||
{ flag cerr;
|
||||
ftnint cunit;
|
||||
char *csta;
|
||||
} cllist;
|
||||
|
||||
/*rewind, backspace, endfile*/
|
||||
typedef struct
|
||||
{ flag aerr;
|
||||
ftnint aunit;
|
||||
} alist;
|
||||
|
||||
/* inquire */
|
||||
typedef struct
|
||||
{ flag inerr;
|
||||
ftnint inunit;
|
||||
char *infile;
|
||||
ftnlen infilen;
|
||||
ftnint *inex; /*parameters in standard's order*/
|
||||
ftnint *inopen;
|
||||
ftnint *innum;
|
||||
ftnint *innamed;
|
||||
char *inname;
|
||||
ftnlen innamlen;
|
||||
char *inacc;
|
||||
ftnlen inacclen;
|
||||
char *inseq;
|
||||
ftnlen inseqlen;
|
||||
char *indir;
|
||||
ftnlen indirlen;
|
||||
char *infmt;
|
||||
ftnlen infmtlen;
|
||||
char *inform;
|
||||
ftnint informlen;
|
||||
char *inunf;
|
||||
ftnlen inunflen;
|
||||
ftnint *inrecl;
|
||||
ftnint *innrec;
|
||||
char *inblank;
|
||||
ftnlen inblanklen;
|
||||
} inlist;
|
||||
|
||||
#define VOID void
|
||||
|
||||
union Multitype { /* for multiple entry points */
|
||||
integer1 g;
|
||||
shortint h;
|
||||
integer i;
|
||||
/* longint j; */
|
||||
real r;
|
||||
doublereal d;
|
||||
complex c;
|
||||
doublecomplex z;
|
||||
};
|
||||
|
||||
typedef union Multitype Multitype;
|
||||
|
||||
struct Vardesc { /* for Namelist */
|
||||
char *name;
|
||||
char *addr;
|
||||
ftnlen *dims;
|
||||
int type;
|
||||
};
|
||||
typedef struct Vardesc Vardesc;
|
||||
|
||||
struct Namelist {
|
||||
char *name;
|
||||
Vardesc **vars;
|
||||
int nvars;
|
||||
};
|
||||
typedef struct Namelist Namelist;
|
||||
|
||||
#define abs(x) ((x) >= 0 ? (x) : -(x))
|
||||
#define dabs(x) (fabs(x))
|
||||
#define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
|
||||
#define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
|
||||
#define dmin(a,b) (f2cmin(a,b))
|
||||
#define dmax(a,b) (f2cmax(a,b))
|
||||
#define bit_test(a,b) ((a) >> (b) & 1)
|
||||
#define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
|
||||
#define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
|
||||
|
||||
#define abort_() { sig_die("Fortran abort routine called", 1); }
|
||||
#define c_abs(z) (cabsf(Cf(z)))
|
||||
#define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
|
||||
#ifdef _MSC_VER
|
||||
#define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
|
||||
#define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/Cd(b)._Val[1]);}
|
||||
#else
|
||||
#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
|
||||
#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
|
||||
#endif
|
||||
#define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
|
||||
#define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
|
||||
#define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
|
||||
//#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
|
||||
#define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
|
||||
#define d_abs(x) (fabs(*(x)))
|
||||
#define d_acos(x) (acos(*(x)))
|
||||
#define d_asin(x) (asin(*(x)))
|
||||
#define d_atan(x) (atan(*(x)))
|
||||
#define d_atn2(x, y) (atan2(*(x),*(y)))
|
||||
#define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
|
||||
#define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
|
||||
#define d_cos(x) (cos(*(x)))
|
||||
#define d_cosh(x) (cosh(*(x)))
|
||||
#define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
|
||||
#define d_exp(x) (exp(*(x)))
|
||||
#define d_imag(z) (cimag(Cd(z)))
|
||||
#define r_imag(z) (cimagf(Cf(z)))
|
||||
#define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define d_log(x) (log(*(x)))
|
||||
#define d_mod(x, y) (fmod(*(x), *(y)))
|
||||
#define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
|
||||
#define d_nint(x) u_nint(*(x))
|
||||
#define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
|
||||
#define d_sign(a,b) u_sign(*(a),*(b))
|
||||
#define r_sign(a,b) u_sign(*(a),*(b))
|
||||
#define d_sin(x) (sin(*(x)))
|
||||
#define d_sinh(x) (sinh(*(x)))
|
||||
#define d_sqrt(x) (sqrt(*(x)))
|
||||
#define d_tan(x) (tan(*(x)))
|
||||
#define d_tanh(x) (tanh(*(x)))
|
||||
#define i_abs(x) abs(*(x))
|
||||
#define i_dnnt(x) ((integer)u_nint(*(x)))
|
||||
#define i_len(s, n) (n)
|
||||
#define i_nint(x) ((integer)u_nint(*(x)))
|
||||
#define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
|
||||
#define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
|
||||
#define pow_si(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_ri(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_di(B,E) dpow_ui(*(B),*(E))
|
||||
#define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
|
||||
#define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
|
||||
#define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
|
||||
#define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
|
||||
#define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
|
||||
#define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
|
||||
#define sig_die(s, kill) { exit(1); }
|
||||
#define s_stop(s, n) {exit(0);}
|
||||
#define z_abs(z) (cabs(Cd(z)))
|
||||
#define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
|
||||
#define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
|
||||
#define myexit_() break;
|
||||
#define mycycle_() continue;
|
||||
#define myceiling_(w) {ceil(w)}
|
||||
#define myhuge_(w) {HUGE_VAL}
|
||||
#define mymaxloc_(w,s,e,n) dmaxloc_(w,*(s),*(e),n)
|
||||
|
||||
/* procedure parameter types for -A and -C++ */
|
||||
|
||||
#define F2C_proc_par_types 1
|
||||
#ifdef __cplusplus
|
||||
typedef logical (*L_fp)(...);
|
||||
#else
|
||||
typedef logical (*L_fp)();
|
||||
#endif
|
||||
|
||||
/* -- translated by f2c (version 20000121).
|
||||
You must link the resulting object file with the libraries:
|
||||
-lf2c -lm (in that order)
|
||||
*/
|
||||
|
||||
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static complex c_b1 = {1.f,0.f};
|
||||
|
||||
/* > \brief \b CGEQRS */
|
||||
|
||||
/* =========== DOCUMENTATION =========== */
|
||||
|
||||
/* Online html documentation available at */
|
||||
/* http://www.netlib.org/lapack/explore-html/ */
|
||||
|
||||
/* Definition: */
|
||||
/* =========== */
|
||||
|
||||
/* SUBROUTINE CGEQRS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, */
|
||||
/* INFO ) */
|
||||
|
||||
/* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS */
|
||||
/* COMPLEX A( LDA, * ), B( LDB, * ), TAU( * ), */
|
||||
/* $ WORK( LWORK ) */
|
||||
|
||||
|
||||
/* > \par Purpose: */
|
||||
/* ============= */
|
||||
/* > */
|
||||
/* > \verbatim */
|
||||
/* > */
|
||||
/* > Solve the least squares problem */
|
||||
/* > f2cmin || A*X - B || */
|
||||
/* > using the QR factorization */
|
||||
/* > A = Q*R */
|
||||
/* > computed by CGEQRF. */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Arguments: */
|
||||
/* ========== */
|
||||
|
||||
/* > \param[in] M */
|
||||
/* > \verbatim */
|
||||
/* > M is INTEGER */
|
||||
/* > The number of rows of the matrix A. M >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] N */
|
||||
/* > \verbatim */
|
||||
/* > N is INTEGER */
|
||||
/* > The number of columns of the matrix A. M >= N >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] NRHS */
|
||||
/* > \verbatim */
|
||||
/* > NRHS is INTEGER */
|
||||
/* > The number of columns of B. NRHS >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] A */
|
||||
/* > \verbatim */
|
||||
/* > A is COMPLEX array, dimension (LDA,N) */
|
||||
/* > Details of the QR factorization of the original matrix A as */
|
||||
/* > returned by CGEQRF. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDA */
|
||||
/* > \verbatim */
|
||||
/* > LDA is INTEGER */
|
||||
/* > The leading dimension of the array A. LDA >= M. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] TAU */
|
||||
/* > \verbatim */
|
||||
/* > TAU is COMPLEX array, dimension (N) */
|
||||
/* > Details of the orthogonal matrix Q. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in,out] B */
|
||||
/* > \verbatim */
|
||||
/* > B is COMPLEX array, dimension (LDB,NRHS) */
|
||||
/* > On entry, the m-by-nrhs right hand side matrix B. */
|
||||
/* > On exit, the n-by-nrhs solution matrix X. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDB */
|
||||
/* > \verbatim */
|
||||
/* > LDB is INTEGER */
|
||||
/* > The leading dimension of the array B. LDB >= M. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] WORK */
|
||||
/* > \verbatim */
|
||||
/* > WORK is COMPLEX array, dimension (LWORK) */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LWORK */
|
||||
/* > \verbatim */
|
||||
/* > LWORK is INTEGER */
|
||||
/* > The length of the array WORK. LWORK must be at least NRHS, */
|
||||
/* > and should be at least NRHS*NB, where NB is the block size */
|
||||
/* > for this environment. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] INFO */
|
||||
/* > \verbatim */
|
||||
/* > INFO is INTEGER */
|
||||
/* > = 0: successful exit */
|
||||
/* > < 0: if INFO = -i, the i-th argument had an illegal value */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Authors: */
|
||||
/* ======== */
|
||||
|
||||
/* > \author Univ. of Tennessee */
|
||||
/* > \author Univ. of California Berkeley */
|
||||
/* > \author Univ. of Colorado Denver */
|
||||
/* > \author NAG Ltd. */
|
||||
|
||||
/* > \ingroup complex_lin */
|
||||
|
||||
/* ===================================================================== */
|
||||
/* Subroutine */ int cgeqrs_(integer *m, integer *n, integer *nrhs, complex *
|
||||
a, integer *lda, complex *tau, complex *b, integer *ldb, complex *
|
||||
work, integer *lwork, integer *info)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, b_dim1, b_offset, i__1;
|
||||
|
||||
/* Local variables */
|
||||
extern /* Subroutine */ int ctrsm_(char *, char *, char *, char *,
|
||||
integer *, integer *, complex *, complex *, integer *, complex *,
|
||||
integer *), xerbla_(char *,
|
||||
integer *), cunmqr_(char *, char *, integer *, integer *,
|
||||
integer *, complex *, integer *, complex *, complex *, integer *,
|
||||
complex *, integer *, integer *);
|
||||
|
||||
|
||||
/* -- LAPACK test routine -- */
|
||||
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
|
||||
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
|
||||
|
||||
|
||||
/* ===================================================================== */
|
||||
|
||||
|
||||
/* Test the input arguments. */
|
||||
|
||||
/* Parameter adjustments */
|
||||
a_dim1 = *lda;
|
||||
a_offset = 1 + a_dim1 * 1;
|
||||
a -= a_offset;
|
||||
--tau;
|
||||
b_dim1 = *ldb;
|
||||
b_offset = 1 + b_dim1 * 1;
|
||||
b -= b_offset;
|
||||
--work;
|
||||
|
||||
/* Function Body */
|
||||
*info = 0;
|
||||
if (*m < 0) {
|
||||
*info = -1;
|
||||
} else if (*n < 0 || *n > *m) {
|
||||
*info = -2;
|
||||
} else if (*nrhs < 0) {
|
||||
*info = -3;
|
||||
} else if (*lda < f2cmax(1,*m)) {
|
||||
*info = -5;
|
||||
} else if (*ldb < f2cmax(1,*m)) {
|
||||
*info = -8;
|
||||
} else if (*lwork < 1 || *lwork < *nrhs && *m > 0 && *n > 0) {
|
||||
*info = -10;
|
||||
}
|
||||
if (*info != 0) {
|
||||
i__1 = -(*info);
|
||||
xerbla_("CGEQRS", &i__1);
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Quick return if possible */
|
||||
|
||||
if (*n == 0 || *nrhs == 0 || *m == 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* B := Q' * B */
|
||||
|
||||
cunmqr_("Left", "Conjugate transpose", m, nrhs, n, &a[a_offset], lda, &
|
||||
tau[1], &b[b_offset], ldb, &work[1], lwork, info);
|
||||
|
||||
/* Solve R*X = B(1:n,:) */
|
||||
|
||||
ctrsm_("Left", "Upper", "No transpose", "Non-unit", n, nrhs, &c_b1, &a[
|
||||
a_offset], lda, &b[b_offset], ldb);
|
||||
|
||||
return 0;
|
||||
|
||||
/* End of CGEQRS */
|
||||
|
||||
} /* cgeqrs_ */
|
||||
|
|
@ -0,0 +1,480 @@
|
|||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <stdio.h>
|
||||
#include <complex.h>
|
||||
#ifdef complex
|
||||
#undef complex
|
||||
#endif
|
||||
#ifdef I
|
||||
#undef I
|
||||
#endif
|
||||
|
||||
#if defined(_WIN64)
|
||||
typedef long long BLASLONG;
|
||||
typedef unsigned long long BLASULONG;
|
||||
#else
|
||||
typedef long BLASLONG;
|
||||
typedef unsigned long BLASULONG;
|
||||
#endif
|
||||
|
||||
#ifdef LAPACK_ILP64
|
||||
typedef BLASLONG blasint;
|
||||
#if defined(_WIN64)
|
||||
#define blasabs(x) llabs(x)
|
||||
#else
|
||||
#define blasabs(x) labs(x)
|
||||
#endif
|
||||
#else
|
||||
typedef int blasint;
|
||||
#define blasabs(x) abs(x)
|
||||
#endif
|
||||
|
||||
typedef blasint integer;
|
||||
|
||||
typedef unsigned int uinteger;
|
||||
typedef char *address;
|
||||
typedef short int shortint;
|
||||
typedef float real;
|
||||
typedef double doublereal;
|
||||
typedef struct { real r, i; } complex;
|
||||
typedef struct { doublereal r, i; } doublecomplex;
|
||||
#ifdef _MSC_VER
|
||||
static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
|
||||
static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
|
||||
static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
|
||||
static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
|
||||
#else
|
||||
static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
|
||||
static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
|
||||
#endif
|
||||
#define pCf(z) (*_pCf(z))
|
||||
#define pCd(z) (*_pCd(z))
|
||||
typedef int logical;
|
||||
typedef short int shortlogical;
|
||||
typedef char logical1;
|
||||
typedef char integer1;
|
||||
|
||||
#define TRUE_ (1)
|
||||
#define FALSE_ (0)
|
||||
|
||||
/* Extern is for use with -E */
|
||||
#ifndef Extern
|
||||
#define Extern extern
|
||||
#endif
|
||||
|
||||
/* I/O stuff */
|
||||
|
||||
typedef int flag;
|
||||
typedef int ftnlen;
|
||||
typedef int ftnint;
|
||||
|
||||
/*external read, write*/
|
||||
typedef struct
|
||||
{ flag cierr;
|
||||
ftnint ciunit;
|
||||
flag ciend;
|
||||
char *cifmt;
|
||||
ftnint cirec;
|
||||
} cilist;
|
||||
|
||||
/*internal read, write*/
|
||||
typedef struct
|
||||
{ flag icierr;
|
||||
char *iciunit;
|
||||
flag iciend;
|
||||
char *icifmt;
|
||||
ftnint icirlen;
|
||||
ftnint icirnum;
|
||||
} icilist;
|
||||
|
||||
/*open*/
|
||||
typedef struct
|
||||
{ flag oerr;
|
||||
ftnint ounit;
|
||||
char *ofnm;
|
||||
ftnlen ofnmlen;
|
||||
char *osta;
|
||||
char *oacc;
|
||||
char *ofm;
|
||||
ftnint orl;
|
||||
char *oblnk;
|
||||
} olist;
|
||||
|
||||
/*close*/
|
||||
typedef struct
|
||||
{ flag cerr;
|
||||
ftnint cunit;
|
||||
char *csta;
|
||||
} cllist;
|
||||
|
||||
/*rewind, backspace, endfile*/
|
||||
typedef struct
|
||||
{ flag aerr;
|
||||
ftnint aunit;
|
||||
} alist;
|
||||
|
||||
/* inquire */
|
||||
typedef struct
|
||||
{ flag inerr;
|
||||
ftnint inunit;
|
||||
char *infile;
|
||||
ftnlen infilen;
|
||||
ftnint *inex; /*parameters in standard's order*/
|
||||
ftnint *inopen;
|
||||
ftnint *innum;
|
||||
ftnint *innamed;
|
||||
char *inname;
|
||||
ftnlen innamlen;
|
||||
char *inacc;
|
||||
ftnlen inacclen;
|
||||
char *inseq;
|
||||
ftnlen inseqlen;
|
||||
char *indir;
|
||||
ftnlen indirlen;
|
||||
char *infmt;
|
||||
ftnlen infmtlen;
|
||||
char *inform;
|
||||
ftnint informlen;
|
||||
char *inunf;
|
||||
ftnlen inunflen;
|
||||
ftnint *inrecl;
|
||||
ftnint *innrec;
|
||||
char *inblank;
|
||||
ftnlen inblanklen;
|
||||
} inlist;
|
||||
|
||||
#define VOID void
|
||||
|
||||
union Multitype { /* for multiple entry points */
|
||||
integer1 g;
|
||||
shortint h;
|
||||
integer i;
|
||||
/* longint j; */
|
||||
real r;
|
||||
doublereal d;
|
||||
complex c;
|
||||
doublecomplex z;
|
||||
};
|
||||
|
||||
typedef union Multitype Multitype;
|
||||
|
||||
struct Vardesc { /* for Namelist */
|
||||
char *name;
|
||||
char *addr;
|
||||
ftnlen *dims;
|
||||
int type;
|
||||
};
|
||||
typedef struct Vardesc Vardesc;
|
||||
|
||||
struct Namelist {
|
||||
char *name;
|
||||
Vardesc **vars;
|
||||
int nvars;
|
||||
};
|
||||
typedef struct Namelist Namelist;
|
||||
|
||||
#define abs(x) ((x) >= 0 ? (x) : -(x))
|
||||
#define dabs(x) (fabs(x))
|
||||
#define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
|
||||
#define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
|
||||
#define dmin(a,b) (f2cmin(a,b))
|
||||
#define dmax(a,b) (f2cmax(a,b))
|
||||
#define bit_test(a,b) ((a) >> (b) & 1)
|
||||
#define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
|
||||
#define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
|
||||
|
||||
#define abort_() { sig_die("Fortran abort routine called", 1); }
|
||||
#define c_abs(z) (cabsf(Cf(z)))
|
||||
#define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
|
||||
#ifdef _MSC_VER
|
||||
#define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
|
||||
#define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/Cd(b)._Val[1]);}
|
||||
#else
|
||||
#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
|
||||
#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
|
||||
#endif
|
||||
#define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
|
||||
#define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
|
||||
#define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
|
||||
//#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
|
||||
#define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
|
||||
#define d_abs(x) (fabs(*(x)))
|
||||
#define d_acos(x) (acos(*(x)))
|
||||
#define d_asin(x) (asin(*(x)))
|
||||
#define d_atan(x) (atan(*(x)))
|
||||
#define d_atn2(x, y) (atan2(*(x),*(y)))
|
||||
#define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
|
||||
#define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
|
||||
#define d_cos(x) (cos(*(x)))
|
||||
#define d_cosh(x) (cosh(*(x)))
|
||||
#define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
|
||||
#define d_exp(x) (exp(*(x)))
|
||||
#define d_imag(z) (cimag(Cd(z)))
|
||||
#define r_imag(z) (cimagf(Cf(z)))
|
||||
#define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define d_log(x) (log(*(x)))
|
||||
#define d_mod(x, y) (fmod(*(x), *(y)))
|
||||
#define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
|
||||
#define d_nint(x) u_nint(*(x))
|
||||
#define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
|
||||
#define d_sign(a,b) u_sign(*(a),*(b))
|
||||
#define r_sign(a,b) u_sign(*(a),*(b))
|
||||
#define d_sin(x) (sin(*(x)))
|
||||
#define d_sinh(x) (sinh(*(x)))
|
||||
#define d_sqrt(x) (sqrt(*(x)))
|
||||
#define d_tan(x) (tan(*(x)))
|
||||
#define d_tanh(x) (tanh(*(x)))
|
||||
#define i_abs(x) abs(*(x))
|
||||
#define i_dnnt(x) ((integer)u_nint(*(x)))
|
||||
#define i_len(s, n) (n)
|
||||
#define i_nint(x) ((integer)u_nint(*(x)))
|
||||
#define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
|
||||
#define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
|
||||
#define pow_si(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_ri(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_di(B,E) dpow_ui(*(B),*(E))
|
||||
#define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
|
||||
#define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
|
||||
#define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
|
||||
#define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
|
||||
#define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
|
||||
#define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
|
||||
#define sig_die(s, kill) { exit(1); }
|
||||
#define s_stop(s, n) {exit(0);}
|
||||
#define z_abs(z) (cabs(Cd(z)))
|
||||
#define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
|
||||
#define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
|
||||
#define myexit_() break;
|
||||
#define mycycle_() continue;
|
||||
#define myceiling_(w) {ceil(w)}
|
||||
#define myhuge_(w) {HUGE_VAL}
|
||||
#define mymaxloc_(w,s,e,n) dmaxloc_(w,*(s),*(e),n)
|
||||
|
||||
/* procedure parameter types for -A and -C++ */
|
||||
|
||||
#define F2C_proc_par_types 1
|
||||
#ifdef __cplusplus
|
||||
typedef logical (*L_fp)(...);
|
||||
#else
|
||||
typedef logical (*L_fp)();
|
||||
#endif
|
||||
|
||||
/* -- translated by f2c (version 20000121).
|
||||
You must link the resulting object file with the libraries:
|
||||
-lf2c -lm (in that order)
|
||||
*/
|
||||
|
||||
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static doublereal c_b7 = 1.;
|
||||
static doublereal c_b9 = 0.;
|
||||
|
||||
/* > \brief \b DGELQS */
|
||||
|
||||
/* =========== DOCUMENTATION =========== */
|
||||
|
||||
/* Online html documentation available at */
|
||||
/* http://www.netlib.org/lapack/explore-html/ */
|
||||
|
||||
/* Definition: */
|
||||
/* =========== */
|
||||
|
||||
/* SUBROUTINE DGELQS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, */
|
||||
/* INFO ) */
|
||||
|
||||
/* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS */
|
||||
/* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), TAU( * ), */
|
||||
/* $ WORK( LWORK ) */
|
||||
|
||||
|
||||
/* > \par Purpose: */
|
||||
/* ============= */
|
||||
/* > */
|
||||
/* > \verbatim */
|
||||
/* > */
|
||||
/* > Compute a minimum-norm solution */
|
||||
/* > f2cmin || A*X - B || */
|
||||
/* > using the LQ factorization */
|
||||
/* > A = L*Q */
|
||||
/* > computed by DGELQF. */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Arguments: */
|
||||
/* ========== */
|
||||
|
||||
/* > \param[in] M */
|
||||
/* > \verbatim */
|
||||
/* > M is INTEGER */
|
||||
/* > The number of rows of the matrix A. M >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] N */
|
||||
/* > \verbatim */
|
||||
/* > N is INTEGER */
|
||||
/* > The number of columns of the matrix A. N >= M >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] NRHS */
|
||||
/* > \verbatim */
|
||||
/* > NRHS is INTEGER */
|
||||
/* > The number of columns of B. NRHS >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] A */
|
||||
/* > \verbatim */
|
||||
/* > A is DOUBLE PRECISION array, dimension (LDA,N) */
|
||||
/* > Details of the LQ factorization of the original matrix A as */
|
||||
/* > returned by DGELQF. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDA */
|
||||
/* > \verbatim */
|
||||
/* > LDA is INTEGER */
|
||||
/* > The leading dimension of the array A. LDA >= M. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] TAU */
|
||||
/* > \verbatim */
|
||||
/* > TAU is DOUBLE PRECISION array, dimension (M) */
|
||||
/* > Details of the orthogonal matrix Q. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in,out] B */
|
||||
/* > \verbatim */
|
||||
/* > B is DOUBLE PRECISION array, dimension (LDB,NRHS) */
|
||||
/* > On entry, the m-by-nrhs right hand side matrix B. */
|
||||
/* > On exit, the n-by-nrhs solution matrix X. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDB */
|
||||
/* > \verbatim */
|
||||
/* > LDB is INTEGER */
|
||||
/* > The leading dimension of the array B. LDB >= N. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] WORK */
|
||||
/* > \verbatim */
|
||||
/* > WORK is DOUBLE PRECISION array, dimension (LWORK) */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LWORK */
|
||||
/* > \verbatim */
|
||||
/* > LWORK is INTEGER */
|
||||
/* > The length of the array WORK. LWORK must be at least NRHS, */
|
||||
/* > and should be at least NRHS*NB, where NB is the block size */
|
||||
/* > for this environment. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] INFO */
|
||||
/* > \verbatim */
|
||||
/* > INFO is INTEGER */
|
||||
/* > = 0: successful exit */
|
||||
/* > < 0: if INFO = -i, the i-th argument had an illegal value */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Authors: */
|
||||
/* ======== */
|
||||
|
||||
/* > \author Univ. of Tennessee */
|
||||
/* > \author Univ. of California Berkeley */
|
||||
/* > \author Univ. of Colorado Denver */
|
||||
/* > \author NAG Ltd. */
|
||||
|
||||
/* > \ingroup double_lin */
|
||||
|
||||
/* ===================================================================== */
|
||||
/* Subroutine */ int dgelqs_(integer *m, integer *n, integer *nrhs,
|
||||
doublereal *a, integer *lda, doublereal *tau, doublereal *b, integer *
|
||||
ldb, doublereal *work, integer *lwork, integer *info)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, b_dim1, b_offset, i__1;
|
||||
|
||||
/* Local variables */
|
||||
extern /* Subroutine */ int dtrsm_(char *, char *, char *, char *,
|
||||
integer *, integer *, doublereal *, doublereal *, integer *,
|
||||
doublereal *, integer *), dlaset_(
|
||||
char *, integer *, integer *, doublereal *, doublereal *,
|
||||
doublereal *, integer *), xerbla_(char *, integer *), dormlq_(char *, char *, integer *, integer *, integer *,
|
||||
doublereal *, integer *, doublereal *, doublereal *, integer *,
|
||||
doublereal *, integer *, integer *);
|
||||
|
||||
|
||||
/* -- LAPACK test routine -- */
|
||||
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
|
||||
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
|
||||
|
||||
|
||||
/* ===================================================================== */
|
||||
|
||||
|
||||
/* Test the input parameters. */
|
||||
|
||||
/* Parameter adjustments */
|
||||
a_dim1 = *lda;
|
||||
a_offset = 1 + a_dim1 * 1;
|
||||
a -= a_offset;
|
||||
--tau;
|
||||
b_dim1 = *ldb;
|
||||
b_offset = 1 + b_dim1 * 1;
|
||||
b -= b_offset;
|
||||
--work;
|
||||
|
||||
/* Function Body */
|
||||
*info = 0;
|
||||
if (*m < 0) {
|
||||
*info = -1;
|
||||
} else if (*n < 0 || *m > *n) {
|
||||
*info = -2;
|
||||
} else if (*nrhs < 0) {
|
||||
*info = -3;
|
||||
} else if (*lda < f2cmax(1,*m)) {
|
||||
*info = -5;
|
||||
} else if (*ldb < f2cmax(1,*n)) {
|
||||
*info = -8;
|
||||
} else if (*lwork < 1 || *lwork < *nrhs && *m > 0 && *n > 0) {
|
||||
*info = -10;
|
||||
}
|
||||
if (*info != 0) {
|
||||
i__1 = -(*info);
|
||||
xerbla_("DGELQS", &i__1);
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Quick return if possible */
|
||||
|
||||
if (*n == 0 || *nrhs == 0 || *m == 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Solve L*X = B(1:m,:) */
|
||||
|
||||
dtrsm_("Left", "Lower", "No transpose", "Non-unit", m, nrhs, &c_b7, &a[
|
||||
a_offset], lda, &b[b_offset], ldb);
|
||||
|
||||
/* Set B(m+1:n,:) to zero */
|
||||
|
||||
if (*m < *n) {
|
||||
i__1 = *n - *m;
|
||||
dlaset_("Full", &i__1, nrhs, &c_b9, &c_b9, &b[*m + 1 + b_dim1], ldb);
|
||||
}
|
||||
|
||||
/* B := Q' * B */
|
||||
|
||||
dormlq_("Left", "Transpose", n, nrhs, m, &a[a_offset], lda, &tau[1], &b[
|
||||
b_offset], ldb, &work[1], lwork, info);
|
||||
|
||||
return 0;
|
||||
|
||||
/* End of DGELQS */
|
||||
|
||||
} /* dgelqs_ */
|
||||
|
|
@ -0,0 +1,471 @@
|
|||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <stdio.h>
|
||||
#include <complex.h>
|
||||
#ifdef complex
|
||||
#undef complex
|
||||
#endif
|
||||
#ifdef I
|
||||
#undef I
|
||||
#endif
|
||||
|
||||
#if defined(_WIN64)
|
||||
typedef long long BLASLONG;
|
||||
typedef unsigned long long BLASULONG;
|
||||
#else
|
||||
typedef long BLASLONG;
|
||||
typedef unsigned long BLASULONG;
|
||||
#endif
|
||||
|
||||
#ifdef LAPACK_ILP64
|
||||
typedef BLASLONG blasint;
|
||||
#if defined(_WIN64)
|
||||
#define blasabs(x) llabs(x)
|
||||
#else
|
||||
#define blasabs(x) labs(x)
|
||||
#endif
|
||||
#else
|
||||
typedef int blasint;
|
||||
#define blasabs(x) abs(x)
|
||||
#endif
|
||||
|
||||
typedef blasint integer;
|
||||
|
||||
typedef unsigned int uinteger;
|
||||
typedef char *address;
|
||||
typedef short int shortint;
|
||||
typedef float real;
|
||||
typedef double doublereal;
|
||||
typedef struct { real r, i; } complex;
|
||||
typedef struct { doublereal r, i; } doublecomplex;
|
||||
#ifdef _MSC_VER
|
||||
static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
|
||||
static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
|
||||
static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
|
||||
static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
|
||||
#else
|
||||
static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
|
||||
static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
|
||||
#endif
|
||||
#define pCf(z) (*_pCf(z))
|
||||
#define pCd(z) (*_pCd(z))
|
||||
typedef int logical;
|
||||
typedef short int shortlogical;
|
||||
typedef char logical1;
|
||||
typedef char integer1;
|
||||
|
||||
#define TRUE_ (1)
|
||||
#define FALSE_ (0)
|
||||
|
||||
/* Extern is for use with -E */
|
||||
#ifndef Extern
|
||||
#define Extern extern
|
||||
#endif
|
||||
|
||||
/* I/O stuff */
|
||||
|
||||
typedef int flag;
|
||||
typedef int ftnlen;
|
||||
typedef int ftnint;
|
||||
|
||||
/*external read, write*/
|
||||
typedef struct
|
||||
{ flag cierr;
|
||||
ftnint ciunit;
|
||||
flag ciend;
|
||||
char *cifmt;
|
||||
ftnint cirec;
|
||||
} cilist;
|
||||
|
||||
/*internal read, write*/
|
||||
typedef struct
|
||||
{ flag icierr;
|
||||
char *iciunit;
|
||||
flag iciend;
|
||||
char *icifmt;
|
||||
ftnint icirlen;
|
||||
ftnint icirnum;
|
||||
} icilist;
|
||||
|
||||
/*open*/
|
||||
typedef struct
|
||||
{ flag oerr;
|
||||
ftnint ounit;
|
||||
char *ofnm;
|
||||
ftnlen ofnmlen;
|
||||
char *osta;
|
||||
char *oacc;
|
||||
char *ofm;
|
||||
ftnint orl;
|
||||
char *oblnk;
|
||||
} olist;
|
||||
|
||||
/*close*/
|
||||
typedef struct
|
||||
{ flag cerr;
|
||||
ftnint cunit;
|
||||
char *csta;
|
||||
} cllist;
|
||||
|
||||
/*rewind, backspace, endfile*/
|
||||
typedef struct
|
||||
{ flag aerr;
|
||||
ftnint aunit;
|
||||
} alist;
|
||||
|
||||
/* inquire */
|
||||
typedef struct
|
||||
{ flag inerr;
|
||||
ftnint inunit;
|
||||
char *infile;
|
||||
ftnlen infilen;
|
||||
ftnint *inex; /*parameters in standard's order*/
|
||||
ftnint *inopen;
|
||||
ftnint *innum;
|
||||
ftnint *innamed;
|
||||
char *inname;
|
||||
ftnlen innamlen;
|
||||
char *inacc;
|
||||
ftnlen inacclen;
|
||||
char *inseq;
|
||||
ftnlen inseqlen;
|
||||
char *indir;
|
||||
ftnlen indirlen;
|
||||
char *infmt;
|
||||
ftnlen infmtlen;
|
||||
char *inform;
|
||||
ftnint informlen;
|
||||
char *inunf;
|
||||
ftnlen inunflen;
|
||||
ftnint *inrecl;
|
||||
ftnint *innrec;
|
||||
char *inblank;
|
||||
ftnlen inblanklen;
|
||||
} inlist;
|
||||
|
||||
#define VOID void
|
||||
|
||||
union Multitype { /* for multiple entry points */
|
||||
integer1 g;
|
||||
shortint h;
|
||||
integer i;
|
||||
/* longint j; */
|
||||
real r;
|
||||
doublereal d;
|
||||
complex c;
|
||||
doublecomplex z;
|
||||
};
|
||||
|
||||
typedef union Multitype Multitype;
|
||||
|
||||
struct Vardesc { /* for Namelist */
|
||||
char *name;
|
||||
char *addr;
|
||||
ftnlen *dims;
|
||||
int type;
|
||||
};
|
||||
typedef struct Vardesc Vardesc;
|
||||
|
||||
struct Namelist {
|
||||
char *name;
|
||||
Vardesc **vars;
|
||||
int nvars;
|
||||
};
|
||||
typedef struct Namelist Namelist;
|
||||
|
||||
#define abs(x) ((x) >= 0 ? (x) : -(x))
|
||||
#define dabs(x) (fabs(x))
|
||||
#define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
|
||||
#define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
|
||||
#define dmin(a,b) (f2cmin(a,b))
|
||||
#define dmax(a,b) (f2cmax(a,b))
|
||||
#define bit_test(a,b) ((a) >> (b) & 1)
|
||||
#define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
|
||||
#define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
|
||||
|
||||
#define abort_() { sig_die("Fortran abort routine called", 1); }
|
||||
#define c_abs(z) (cabsf(Cf(z)))
|
||||
#define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
|
||||
#ifdef _MSC_VER
|
||||
#define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
|
||||
#define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/Cd(b)._Val[1]);}
|
||||
#else
|
||||
#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
|
||||
#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
|
||||
#endif
|
||||
#define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
|
||||
#define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
|
||||
#define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
|
||||
//#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
|
||||
#define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
|
||||
#define d_abs(x) (fabs(*(x)))
|
||||
#define d_acos(x) (acos(*(x)))
|
||||
#define d_asin(x) (asin(*(x)))
|
||||
#define d_atan(x) (atan(*(x)))
|
||||
#define d_atn2(x, y) (atan2(*(x),*(y)))
|
||||
#define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
|
||||
#define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
|
||||
#define d_cos(x) (cos(*(x)))
|
||||
#define d_cosh(x) (cosh(*(x)))
|
||||
#define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
|
||||
#define d_exp(x) (exp(*(x)))
|
||||
#define d_imag(z) (cimag(Cd(z)))
|
||||
#define r_imag(z) (cimagf(Cf(z)))
|
||||
#define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define d_log(x) (log(*(x)))
|
||||
#define d_mod(x, y) (fmod(*(x), *(y)))
|
||||
#define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
|
||||
#define d_nint(x) u_nint(*(x))
|
||||
#define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
|
||||
#define d_sign(a,b) u_sign(*(a),*(b))
|
||||
#define r_sign(a,b) u_sign(*(a),*(b))
|
||||
#define d_sin(x) (sin(*(x)))
|
||||
#define d_sinh(x) (sinh(*(x)))
|
||||
#define d_sqrt(x) (sqrt(*(x)))
|
||||
#define d_tan(x) (tan(*(x)))
|
||||
#define d_tanh(x) (tanh(*(x)))
|
||||
#define i_abs(x) abs(*(x))
|
||||
#define i_dnnt(x) ((integer)u_nint(*(x)))
|
||||
#define i_len(s, n) (n)
|
||||
#define i_nint(x) ((integer)u_nint(*(x)))
|
||||
#define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
|
||||
#define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
|
||||
#define pow_si(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_ri(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_di(B,E) dpow_ui(*(B),*(E))
|
||||
#define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
|
||||
#define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
|
||||
#define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
|
||||
#define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
|
||||
#define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
|
||||
#define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
|
||||
#define sig_die(s, kill) { exit(1); }
|
||||
#define s_stop(s, n) {exit(0);}
|
||||
#define z_abs(z) (cabs(Cd(z)))
|
||||
#define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
|
||||
#define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
|
||||
#define myexit_() break;
|
||||
#define mycycle_() continue;
|
||||
#define myceiling_(w) {ceil(w)}
|
||||
#define myhuge_(w) {HUGE_VAL}
|
||||
#define mymaxloc_(w,s,e,n) dmaxloc_(w,*(s),*(e),n)
|
||||
|
||||
/* procedure parameter types for -A and -C++ */
|
||||
|
||||
#define F2C_proc_par_types 1
|
||||
#ifdef __cplusplus
|
||||
typedef logical (*L_fp)(...);
|
||||
#else
|
||||
typedef logical (*L_fp)();
|
||||
#endif
|
||||
|
||||
/* -- translated by f2c (version 20000121).
|
||||
You must link the resulting object file with the libraries:
|
||||
-lf2c -lm (in that order)
|
||||
*/
|
||||
|
||||
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static doublereal c_b9 = 1.;
|
||||
|
||||
/* > \brief \b DGEQRS */
|
||||
|
||||
/* =========== DOCUMENTATION =========== */
|
||||
|
||||
/* Online html documentation available at */
|
||||
/* http://www.netlib.org/lapack/explore-html/ */
|
||||
|
||||
/* Definition: */
|
||||
/* =========== */
|
||||
|
||||
/* SUBROUTINE DGEQRS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, */
|
||||
/* INFO ) */
|
||||
|
||||
/* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS */
|
||||
/* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), TAU( * ), */
|
||||
/* $ WORK( LWORK ) */
|
||||
|
||||
|
||||
/* > \par Purpose: */
|
||||
/* ============= */
|
||||
/* > */
|
||||
/* > \verbatim */
|
||||
/* > */
|
||||
/* > Solve the least squares problem */
|
||||
/* > f2cmin || A*X - B || */
|
||||
/* > using the QR factorization */
|
||||
/* > A = Q*R */
|
||||
/* > computed by DGEQRF. */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Arguments: */
|
||||
/* ========== */
|
||||
|
||||
/* > \param[in] M */
|
||||
/* > \verbatim */
|
||||
/* > M is INTEGER */
|
||||
/* > The number of rows of the matrix A. M >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] N */
|
||||
/* > \verbatim */
|
||||
/* > N is INTEGER */
|
||||
/* > The number of columns of the matrix A. M >= N >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] NRHS */
|
||||
/* > \verbatim */
|
||||
/* > NRHS is INTEGER */
|
||||
/* > The number of columns of B. NRHS >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] A */
|
||||
/* > \verbatim */
|
||||
/* > A is DOUBLE PRECISION array, dimension (LDA,N) */
|
||||
/* > Details of the QR factorization of the original matrix A as */
|
||||
/* > returned by DGEQRF. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDA */
|
||||
/* > \verbatim */
|
||||
/* > LDA is INTEGER */
|
||||
/* > The leading dimension of the array A. LDA >= M. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] TAU */
|
||||
/* > \verbatim */
|
||||
/* > TAU is DOUBLE PRECISION array, dimension (N) */
|
||||
/* > Details of the orthogonal matrix Q. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in,out] B */
|
||||
/* > \verbatim */
|
||||
/* > B is DOUBLE PRECISION array, dimension (LDB,NRHS) */
|
||||
/* > On entry, the m-by-nrhs right hand side matrix B. */
|
||||
/* > On exit, the n-by-nrhs solution matrix X. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDB */
|
||||
/* > \verbatim */
|
||||
/* > LDB is INTEGER */
|
||||
/* > The leading dimension of the array B. LDB >= M. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] WORK */
|
||||
/* > \verbatim */
|
||||
/* > WORK is DOUBLE PRECISION array, dimension (LWORK) */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LWORK */
|
||||
/* > \verbatim */
|
||||
/* > LWORK is INTEGER */
|
||||
/* > The length of the array WORK. LWORK must be at least NRHS, */
|
||||
/* > and should be at least NRHS*NB, where NB is the block size */
|
||||
/* > for this environment. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] INFO */
|
||||
/* > \verbatim */
|
||||
/* > INFO is INTEGER */
|
||||
/* > = 0: successful exit */
|
||||
/* > < 0: if INFO = -i, the i-th argument had an illegal value */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Authors: */
|
||||
/* ======== */
|
||||
|
||||
/* > \author Univ. of Tennessee */
|
||||
/* > \author Univ. of California Berkeley */
|
||||
/* > \author Univ. of Colorado Denver */
|
||||
/* > \author NAG Ltd. */
|
||||
|
||||
/* > \ingroup double_lin */
|
||||
|
||||
/* ===================================================================== */
|
||||
/* Subroutine */ int dgeqrs_(integer *m, integer *n, integer *nrhs,
|
||||
doublereal *a, integer *lda, doublereal *tau, doublereal *b, integer *
|
||||
ldb, doublereal *work, integer *lwork, integer *info)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, b_dim1, b_offset, i__1;
|
||||
|
||||
/* Local variables */
|
||||
extern /* Subroutine */ int dtrsm_(char *, char *, char *, char *,
|
||||
integer *, integer *, doublereal *, doublereal *, integer *,
|
||||
doublereal *, integer *), xerbla_(
|
||||
char *, integer *), dormqr_(char *, char *, integer *,
|
||||
integer *, integer *, doublereal *, integer *, doublereal *,
|
||||
doublereal *, integer *, doublereal *, integer *, integer *);
|
||||
|
||||
|
||||
/* -- LAPACK test routine -- */
|
||||
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
|
||||
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
|
||||
|
||||
|
||||
/* ===================================================================== */
|
||||
|
||||
|
||||
/* Test the input arguments. */
|
||||
|
||||
/* Parameter adjustments */
|
||||
a_dim1 = *lda;
|
||||
a_offset = 1 + a_dim1 * 1;
|
||||
a -= a_offset;
|
||||
--tau;
|
||||
b_dim1 = *ldb;
|
||||
b_offset = 1 + b_dim1 * 1;
|
||||
b -= b_offset;
|
||||
--work;
|
||||
|
||||
/* Function Body */
|
||||
*info = 0;
|
||||
if (*m < 0) {
|
||||
*info = -1;
|
||||
} else if (*n < 0 || *n > *m) {
|
||||
*info = -2;
|
||||
} else if (*nrhs < 0) {
|
||||
*info = -3;
|
||||
} else if (*lda < f2cmax(1,*m)) {
|
||||
*info = -5;
|
||||
} else if (*ldb < f2cmax(1,*m)) {
|
||||
*info = -8;
|
||||
} else if (*lwork < 1 || *lwork < *nrhs && *m > 0 && *n > 0) {
|
||||
*info = -10;
|
||||
}
|
||||
if (*info != 0) {
|
||||
i__1 = -(*info);
|
||||
xerbla_("DGEQRS", &i__1);
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Quick return if possible */
|
||||
|
||||
if (*n == 0 || *nrhs == 0 || *m == 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* B := Q' * B */
|
||||
|
||||
dormqr_("Left", "Transpose", m, nrhs, n, &a[a_offset], lda, &tau[1], &b[
|
||||
b_offset], ldb, &work[1], lwork, info);
|
||||
|
||||
/* Solve R*X = B(1:n,:) */
|
||||
|
||||
dtrsm_("Left", "Upper", "No transpose", "Non-unit", n, nrhs, &c_b9, &a[
|
||||
a_offset], lda, &b[b_offset], ldb);
|
||||
|
||||
return 0;
|
||||
|
||||
/* End of DGEQRS */
|
||||
|
||||
} /* dgeqrs_ */
|
||||
|
|
@ -0,0 +1,472 @@
|
|||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <stdio.h>
|
||||
#include <complex.h>
|
||||
#ifdef complex
|
||||
#undef complex
|
||||
#endif
|
||||
#ifdef I
|
||||
#undef I
|
||||
#endif
|
||||
|
||||
#if defined(_WIN64)
|
||||
typedef long long BLASLONG;
|
||||
typedef unsigned long long BLASULONG;
|
||||
#else
|
||||
typedef long BLASLONG;
|
||||
typedef unsigned long BLASULONG;
|
||||
#endif
|
||||
|
||||
#ifdef LAPACK_ILP64
|
||||
typedef BLASLONG blasint;
|
||||
#if defined(_WIN64)
|
||||
#define blasabs(x) llabs(x)
|
||||
#else
|
||||
#define blasabs(x) labs(x)
|
||||
#endif
|
||||
#else
|
||||
typedef int blasint;
|
||||
#define blasabs(x) abs(x)
|
||||
#endif
|
||||
|
||||
typedef blasint integer;
|
||||
|
||||
typedef unsigned int uinteger;
|
||||
typedef char *address;
|
||||
typedef short int shortint;
|
||||
typedef float real;
|
||||
typedef double doublereal;
|
||||
typedef struct { real r, i; } complex;
|
||||
typedef struct { doublereal r, i; } doublecomplex;
|
||||
#ifdef _MSC_VER
|
||||
static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
|
||||
static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
|
||||
static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
|
||||
static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
|
||||
#else
|
||||
static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
|
||||
static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
|
||||
#endif
|
||||
#define pCf(z) (*_pCf(z))
|
||||
#define pCd(z) (*_pCd(z))
|
||||
typedef int logical;
|
||||
typedef short int shortlogical;
|
||||
typedef char logical1;
|
||||
typedef char integer1;
|
||||
|
||||
#define TRUE_ (1)
|
||||
#define FALSE_ (0)
|
||||
|
||||
/* Extern is for use with -E */
|
||||
#ifndef Extern
|
||||
#define Extern extern
|
||||
#endif
|
||||
|
||||
/* I/O stuff */
|
||||
|
||||
typedef int flag;
|
||||
typedef int ftnlen;
|
||||
typedef int ftnint;
|
||||
|
||||
/*external read, write*/
|
||||
typedef struct
|
||||
{ flag cierr;
|
||||
ftnint ciunit;
|
||||
flag ciend;
|
||||
char *cifmt;
|
||||
ftnint cirec;
|
||||
} cilist;
|
||||
|
||||
/*internal read, write*/
|
||||
typedef struct
|
||||
{ flag icierr;
|
||||
char *iciunit;
|
||||
flag iciend;
|
||||
char *icifmt;
|
||||
ftnint icirlen;
|
||||
ftnint icirnum;
|
||||
} icilist;
|
||||
|
||||
/*open*/
|
||||
typedef struct
|
||||
{ flag oerr;
|
||||
ftnint ounit;
|
||||
char *ofnm;
|
||||
ftnlen ofnmlen;
|
||||
char *osta;
|
||||
char *oacc;
|
||||
char *ofm;
|
||||
ftnint orl;
|
||||
char *oblnk;
|
||||
} olist;
|
||||
|
||||
/*close*/
|
||||
typedef struct
|
||||
{ flag cerr;
|
||||
ftnint cunit;
|
||||
char *csta;
|
||||
} cllist;
|
||||
|
||||
/*rewind, backspace, endfile*/
|
||||
typedef struct
|
||||
{ flag aerr;
|
||||
ftnint aunit;
|
||||
} alist;
|
||||
|
||||
/* inquire */
|
||||
typedef struct
|
||||
{ flag inerr;
|
||||
ftnint inunit;
|
||||
char *infile;
|
||||
ftnlen infilen;
|
||||
ftnint *inex; /*parameters in standard's order*/
|
||||
ftnint *inopen;
|
||||
ftnint *innum;
|
||||
ftnint *innamed;
|
||||
char *inname;
|
||||
ftnlen innamlen;
|
||||
char *inacc;
|
||||
ftnlen inacclen;
|
||||
char *inseq;
|
||||
ftnlen inseqlen;
|
||||
char *indir;
|
||||
ftnlen indirlen;
|
||||
char *infmt;
|
||||
ftnlen infmtlen;
|
||||
char *inform;
|
||||
ftnint informlen;
|
||||
char *inunf;
|
||||
ftnlen inunflen;
|
||||
ftnint *inrecl;
|
||||
ftnint *innrec;
|
||||
char *inblank;
|
||||
ftnlen inblanklen;
|
||||
} inlist;
|
||||
|
||||
#define VOID void
|
||||
|
||||
union Multitype { /* for multiple entry points */
|
||||
integer1 g;
|
||||
shortint h;
|
||||
integer i;
|
||||
/* longint j; */
|
||||
real r;
|
||||
doublereal d;
|
||||
complex c;
|
||||
doublecomplex z;
|
||||
};
|
||||
|
||||
typedef union Multitype Multitype;
|
||||
|
||||
struct Vardesc { /* for Namelist */
|
||||
char *name;
|
||||
char *addr;
|
||||
ftnlen *dims;
|
||||
int type;
|
||||
};
|
||||
typedef struct Vardesc Vardesc;
|
||||
|
||||
struct Namelist {
|
||||
char *name;
|
||||
Vardesc **vars;
|
||||
int nvars;
|
||||
};
|
||||
typedef struct Namelist Namelist;
|
||||
|
||||
#define abs(x) ((x) >= 0 ? (x) : -(x))
|
||||
#define dabs(x) (fabs(x))
|
||||
#define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
|
||||
#define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
|
||||
#define dmin(a,b) (f2cmin(a,b))
|
||||
#define dmax(a,b) (f2cmax(a,b))
|
||||
#define bit_test(a,b) ((a) >> (b) & 1)
|
||||
#define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
|
||||
#define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
|
||||
|
||||
#define abort_() { sig_die("Fortran abort routine called", 1); }
|
||||
#define c_abs(z) (cabsf(Cf(z)))
|
||||
#define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
|
||||
#ifdef _MSC_VER
|
||||
#define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
|
||||
#define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/Cd(b)._Val[1]);}
|
||||
#else
|
||||
#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
|
||||
#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
|
||||
#endif
|
||||
#define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
|
||||
#define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
|
||||
#define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
|
||||
//#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
|
||||
#define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
|
||||
#define d_abs(x) (fabs(*(x)))
|
||||
#define d_acos(x) (acos(*(x)))
|
||||
#define d_asin(x) (asin(*(x)))
|
||||
#define d_atan(x) (atan(*(x)))
|
||||
#define d_atn2(x, y) (atan2(*(x),*(y)))
|
||||
#define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
|
||||
#define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
|
||||
#define d_cos(x) (cos(*(x)))
|
||||
#define d_cosh(x) (cosh(*(x)))
|
||||
#define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
|
||||
#define d_exp(x) (exp(*(x)))
|
||||
#define d_imag(z) (cimag(Cd(z)))
|
||||
#define r_imag(z) (cimagf(Cf(z)))
|
||||
#define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define d_log(x) (log(*(x)))
|
||||
#define d_mod(x, y) (fmod(*(x), *(y)))
|
||||
#define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
|
||||
#define d_nint(x) u_nint(*(x))
|
||||
#define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
|
||||
#define d_sign(a,b) u_sign(*(a),*(b))
|
||||
#define r_sign(a,b) u_sign(*(a),*(b))
|
||||
#define d_sin(x) (sin(*(x)))
|
||||
#define d_sinh(x) (sinh(*(x)))
|
||||
#define d_sqrt(x) (sqrt(*(x)))
|
||||
#define d_tan(x) (tan(*(x)))
|
||||
#define d_tanh(x) (tanh(*(x)))
|
||||
#define i_abs(x) abs(*(x))
|
||||
#define i_dnnt(x) ((integer)u_nint(*(x)))
|
||||
#define i_len(s, n) (n)
|
||||
#define i_nint(x) ((integer)u_nint(*(x)))
|
||||
#define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
|
||||
#define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
|
||||
#define pow_si(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_ri(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_di(B,E) dpow_ui(*(B),*(E))
|
||||
#define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
|
||||
#define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
|
||||
#define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
|
||||
#define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
|
||||
#define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
|
||||
#define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
|
||||
#define sig_die(s, kill) { exit(1); }
|
||||
#define s_stop(s, n) {exit(0);}
|
||||
#define z_abs(z) (cabs(Cd(z)))
|
||||
#define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
|
||||
#define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
|
||||
#define myexit_() break;
|
||||
#define mycycle_() continue;
|
||||
#define myceiling_(w) {ceil(w)}
|
||||
#define myhuge_(w) {HUGE_VAL}
|
||||
#define mymaxloc_(w,s,e,n) dmaxloc_(w,*(s),*(e),n)
|
||||
|
||||
/* procedure parameter types for -A and -C++ */
|
||||
|
||||
#define F2C_proc_par_types 1
|
||||
#ifdef __cplusplus
|
||||
typedef logical (*L_fp)(...);
|
||||
#else
|
||||
typedef logical (*L_fp)();
|
||||
#endif
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static real c_b7 = 1.f;
|
||||
static real c_b9 = 0.f;
|
||||
|
||||
/* > \brief \b SGELQS */
|
||||
|
||||
/* =========== DOCUMENTATION =========== */
|
||||
|
||||
/* Online html documentation available at */
|
||||
/* http://www.netlib.org/lapack/explore-html/ */
|
||||
|
||||
/* Definition: */
|
||||
/* =========== */
|
||||
|
||||
/* SUBROUTINE SGELQS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, */
|
||||
/* INFO ) */
|
||||
|
||||
/* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS */
|
||||
/* REAL A( LDA, * ), B( LDB, * ), TAU( * ), */
|
||||
/* $ WORK( LWORK ) */
|
||||
|
||||
|
||||
/* > \par Purpose: */
|
||||
/* ============= */
|
||||
/* > */
|
||||
/* > \verbatim */
|
||||
/* > */
|
||||
/* > Compute a minimum-norm solution */
|
||||
/* > f2cmin || A*X - B || */
|
||||
/* > using the LQ factorization */
|
||||
/* > A = L*Q */
|
||||
/* > computed by SGELQF. */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Arguments: */
|
||||
/* ========== */
|
||||
|
||||
/* > \param[in] M */
|
||||
/* > \verbatim */
|
||||
/* > M is INTEGER */
|
||||
/* > The number of rows of the matrix A. M >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] N */
|
||||
/* > \verbatim */
|
||||
/* > N is INTEGER */
|
||||
/* > The number of columns of the matrix A. N >= M >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] NRHS */
|
||||
/* > \verbatim */
|
||||
/* > NRHS is INTEGER */
|
||||
/* > The number of columns of B. NRHS >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] A */
|
||||
/* > \verbatim */
|
||||
/* > A is REAL array, dimension (LDA,N) */
|
||||
/* > Details of the LQ factorization of the original matrix A as */
|
||||
/* > returned by SGELQF. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDA */
|
||||
/* > \verbatim */
|
||||
/* > LDA is INTEGER */
|
||||
/* > The leading dimension of the array A. LDA >= M. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] TAU */
|
||||
/* > \verbatim */
|
||||
/* > TAU is REAL array, dimension (M) */
|
||||
/* > Details of the orthogonal matrix Q. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in,out] B */
|
||||
/* > \verbatim */
|
||||
/* > B is REAL array, dimension (LDB,NRHS) */
|
||||
/* > On entry, the m-by-nrhs right hand side matrix B. */
|
||||
/* > On exit, the n-by-nrhs solution matrix X. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDB */
|
||||
/* > \verbatim */
|
||||
/* > LDB is INTEGER */
|
||||
/* > The leading dimension of the array B. LDB >= N. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] WORK */
|
||||
/* > \verbatim */
|
||||
/* > WORK is REAL array, dimension (LWORK) */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LWORK */
|
||||
/* > \verbatim */
|
||||
/* > LWORK is INTEGER */
|
||||
/* > The length of the array WORK. LWORK must be at least NRHS, */
|
||||
/* > and should be at least NRHS*NB, where NB is the block size */
|
||||
/* > for this environment. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] INFO */
|
||||
/* > \verbatim */
|
||||
/* > INFO is INTEGER */
|
||||
/* > = 0: successful exit */
|
||||
/* > < 0: if INFO = -i, the i-th argument had an illegal value */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Authors: */
|
||||
/* ======== */
|
||||
|
||||
/* > \author Univ. of Tennessee */
|
||||
/* > \author Univ. of California Berkeley */
|
||||
/* > \author Univ. of Colorado Denver */
|
||||
/* > \author NAG Ltd. */
|
||||
|
||||
/* > \ingroup single_lin */
|
||||
|
||||
/* ===================================================================== */
|
||||
/* Subroutine */ int sgelqs_(integer *m, integer *n, integer *nrhs, real *a,
|
||||
integer *lda, real *tau, real *b, integer *ldb, real *work, integer *
|
||||
lwork, integer *info)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, b_dim1, b_offset, i__1;
|
||||
|
||||
/* Local variables */
|
||||
extern /* Subroutine */ int strsm_(char *, char *, char *, char *,
|
||||
integer *, integer *, real *, real *, integer *, real *, integer *
|
||||
), xerbla_(char *, integer *), slaset_(char *, integer *, integer *, real *, real *,
|
||||
real *, integer *), sormlq_(char *, char *, integer *,
|
||||
integer *, integer *, real *, integer *, real *, real *, integer *
|
||||
, real *, integer *, integer *);
|
||||
|
||||
|
||||
/* -- LAPACK test routine -- */
|
||||
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
|
||||
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
|
||||
|
||||
|
||||
/* ===================================================================== */
|
||||
|
||||
|
||||
/* Test the input parameters. */
|
||||
|
||||
/* Parameter adjustments */
|
||||
a_dim1 = *lda;
|
||||
a_offset = 1 + a_dim1 * 1;
|
||||
a -= a_offset;
|
||||
--tau;
|
||||
b_dim1 = *ldb;
|
||||
b_offset = 1 + b_dim1 * 1;
|
||||
b -= b_offset;
|
||||
--work;
|
||||
|
||||
/* Function Body */
|
||||
*info = 0;
|
||||
if (*m < 0) {
|
||||
*info = -1;
|
||||
} else if (*n < 0 || *m > *n) {
|
||||
*info = -2;
|
||||
} else if (*nrhs < 0) {
|
||||
*info = -3;
|
||||
} else if (*lda < f2cmax(1,*m)) {
|
||||
*info = -5;
|
||||
} else if (*ldb < f2cmax(1,*n)) {
|
||||
*info = -8;
|
||||
} else if (*lwork < 1 || *lwork < *nrhs && *m > 0 && *n > 0) {
|
||||
*info = -10;
|
||||
}
|
||||
if (*info != 0) {
|
||||
i__1 = -(*info);
|
||||
xerbla_("SGELQS", &i__1);
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Quick return if possible */
|
||||
|
||||
if (*n == 0 || *nrhs == 0 || *m == 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Solve L*X = B(1:m,:) */
|
||||
|
||||
strsm_("Left", "Lower", "No transpose", "Non-unit", m, nrhs, &c_b7, &a[
|
||||
a_offset], lda, &b[b_offset], ldb);
|
||||
|
||||
/* Set B(m+1:n,:) to zero */
|
||||
|
||||
if (*m < *n) {
|
||||
i__1 = *n - *m;
|
||||
slaset_("Full", &i__1, nrhs, &c_b9, &c_b9, &b[*m + 1 + b_dim1], ldb);
|
||||
}
|
||||
|
||||
/* B := Q' * B */
|
||||
|
||||
sormlq_("Left", "Transpose", n, nrhs, m, &a[a_offset], lda, &tau[1], &b[
|
||||
b_offset], ldb, &work[1], lwork, info);
|
||||
|
||||
return 0;
|
||||
|
||||
/* End of SGELQS */
|
||||
|
||||
} /* sgelqs_ */
|
||||
|
|
@ -0,0 +1,470 @@
|
|||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <stdio.h>
|
||||
#include <complex.h>
|
||||
#ifdef complex
|
||||
#undef complex
|
||||
#endif
|
||||
#ifdef I
|
||||
#undef I
|
||||
#endif
|
||||
|
||||
#if defined(_WIN64)
|
||||
typedef long long BLASLONG;
|
||||
typedef unsigned long long BLASULONG;
|
||||
#else
|
||||
typedef long BLASLONG;
|
||||
typedef unsigned long BLASULONG;
|
||||
#endif
|
||||
|
||||
#ifdef LAPACK_ILP64
|
||||
typedef BLASLONG blasint;
|
||||
#if defined(_WIN64)
|
||||
#define blasabs(x) llabs(x)
|
||||
#else
|
||||
#define blasabs(x) labs(x)
|
||||
#endif
|
||||
#else
|
||||
typedef int blasint;
|
||||
#define blasabs(x) abs(x)
|
||||
#endif
|
||||
|
||||
typedef blasint integer;
|
||||
|
||||
typedef unsigned int uinteger;
|
||||
typedef char *address;
|
||||
typedef short int shortint;
|
||||
typedef float real;
|
||||
typedef double doublereal;
|
||||
typedef struct { real r, i; } complex;
|
||||
typedef struct { doublereal r, i; } doublecomplex;
|
||||
#ifdef _MSC_VER
|
||||
static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
|
||||
static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
|
||||
static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
|
||||
static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
|
||||
#else
|
||||
static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
|
||||
static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
|
||||
#endif
|
||||
#define pCf(z) (*_pCf(z))
|
||||
#define pCd(z) (*_pCd(z))
|
||||
typedef int logical;
|
||||
typedef short int shortlogical;
|
||||
typedef char logical1;
|
||||
typedef char integer1;
|
||||
|
||||
#define TRUE_ (1)
|
||||
#define FALSE_ (0)
|
||||
|
||||
/* Extern is for use with -E */
|
||||
#ifndef Extern
|
||||
#define Extern extern
|
||||
#endif
|
||||
|
||||
/* I/O stuff */
|
||||
|
||||
typedef int flag;
|
||||
typedef int ftnlen;
|
||||
typedef int ftnint;
|
||||
|
||||
/*external read, write*/
|
||||
typedef struct
|
||||
{ flag cierr;
|
||||
ftnint ciunit;
|
||||
flag ciend;
|
||||
char *cifmt;
|
||||
ftnint cirec;
|
||||
} cilist;
|
||||
|
||||
/*internal read, write*/
|
||||
typedef struct
|
||||
{ flag icierr;
|
||||
char *iciunit;
|
||||
flag iciend;
|
||||
char *icifmt;
|
||||
ftnint icirlen;
|
||||
ftnint icirnum;
|
||||
} icilist;
|
||||
|
||||
/*open*/
|
||||
typedef struct
|
||||
{ flag oerr;
|
||||
ftnint ounit;
|
||||
char *ofnm;
|
||||
ftnlen ofnmlen;
|
||||
char *osta;
|
||||
char *oacc;
|
||||
char *ofm;
|
||||
ftnint orl;
|
||||
char *oblnk;
|
||||
} olist;
|
||||
|
||||
/*close*/
|
||||
typedef struct
|
||||
{ flag cerr;
|
||||
ftnint cunit;
|
||||
char *csta;
|
||||
} cllist;
|
||||
|
||||
/*rewind, backspace, endfile*/
|
||||
typedef struct
|
||||
{ flag aerr;
|
||||
ftnint aunit;
|
||||
} alist;
|
||||
|
||||
/* inquire */
|
||||
typedef struct
|
||||
{ flag inerr;
|
||||
ftnint inunit;
|
||||
char *infile;
|
||||
ftnlen infilen;
|
||||
ftnint *inex; /*parameters in standard's order*/
|
||||
ftnint *inopen;
|
||||
ftnint *innum;
|
||||
ftnint *innamed;
|
||||
char *inname;
|
||||
ftnlen innamlen;
|
||||
char *inacc;
|
||||
ftnlen inacclen;
|
||||
char *inseq;
|
||||
ftnlen inseqlen;
|
||||
char *indir;
|
||||
ftnlen indirlen;
|
||||
char *infmt;
|
||||
ftnlen infmtlen;
|
||||
char *inform;
|
||||
ftnint informlen;
|
||||
char *inunf;
|
||||
ftnlen inunflen;
|
||||
ftnint *inrecl;
|
||||
ftnint *innrec;
|
||||
char *inblank;
|
||||
ftnlen inblanklen;
|
||||
} inlist;
|
||||
|
||||
#define VOID void
|
||||
|
||||
union Multitype { /* for multiple entry points */
|
||||
integer1 g;
|
||||
shortint h;
|
||||
integer i;
|
||||
/* longint j; */
|
||||
real r;
|
||||
doublereal d;
|
||||
complex c;
|
||||
doublecomplex z;
|
||||
};
|
||||
|
||||
typedef union Multitype Multitype;
|
||||
|
||||
struct Vardesc { /* for Namelist */
|
||||
char *name;
|
||||
char *addr;
|
||||
ftnlen *dims;
|
||||
int type;
|
||||
};
|
||||
typedef struct Vardesc Vardesc;
|
||||
|
||||
struct Namelist {
|
||||
char *name;
|
||||
Vardesc **vars;
|
||||
int nvars;
|
||||
};
|
||||
typedef struct Namelist Namelist;
|
||||
|
||||
#define abs(x) ((x) >= 0 ? (x) : -(x))
|
||||
#define dabs(x) (fabs(x))
|
||||
#define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
|
||||
#define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
|
||||
#define dmin(a,b) (f2cmin(a,b))
|
||||
#define dmax(a,b) (f2cmax(a,b))
|
||||
#define bit_test(a,b) ((a) >> (b) & 1)
|
||||
#define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
|
||||
#define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
|
||||
|
||||
#define abort_() { sig_die("Fortran abort routine called", 1); }
|
||||
#define c_abs(z) (cabsf(Cf(z)))
|
||||
#define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
|
||||
#ifdef _MSC_VER
|
||||
#define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
|
||||
#define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/Cd(b)._Val[1]);}
|
||||
#else
|
||||
#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
|
||||
#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
|
||||
#endif
|
||||
#define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
|
||||
#define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
|
||||
#define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
|
||||
//#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
|
||||
#define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
|
||||
#define d_abs(x) (fabs(*(x)))
|
||||
#define d_acos(x) (acos(*(x)))
|
||||
#define d_asin(x) (asin(*(x)))
|
||||
#define d_atan(x) (atan(*(x)))
|
||||
#define d_atn2(x, y) (atan2(*(x),*(y)))
|
||||
#define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
|
||||
#define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
|
||||
#define d_cos(x) (cos(*(x)))
|
||||
#define d_cosh(x) (cosh(*(x)))
|
||||
#define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
|
||||
#define d_exp(x) (exp(*(x)))
|
||||
#define d_imag(z) (cimag(Cd(z)))
|
||||
#define r_imag(z) (cimagf(Cf(z)))
|
||||
#define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define d_log(x) (log(*(x)))
|
||||
#define d_mod(x, y) (fmod(*(x), *(y)))
|
||||
#define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
|
||||
#define d_nint(x) u_nint(*(x))
|
||||
#define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
|
||||
#define d_sign(a,b) u_sign(*(a),*(b))
|
||||
#define r_sign(a,b) u_sign(*(a),*(b))
|
||||
#define d_sin(x) (sin(*(x)))
|
||||
#define d_sinh(x) (sinh(*(x)))
|
||||
#define d_sqrt(x) (sqrt(*(x)))
|
||||
#define d_tan(x) (tan(*(x)))
|
||||
#define d_tanh(x) (tanh(*(x)))
|
||||
#define i_abs(x) abs(*(x))
|
||||
#define i_dnnt(x) ((integer)u_nint(*(x)))
|
||||
#define i_len(s, n) (n)
|
||||
#define i_nint(x) ((integer)u_nint(*(x)))
|
||||
#define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
|
||||
#define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
|
||||
#define pow_si(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_ri(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_di(B,E) dpow_ui(*(B),*(E))
|
||||
#define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
|
||||
#define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
|
||||
#define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
|
||||
#define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
|
||||
#define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
|
||||
#define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
|
||||
#define sig_die(s, kill) { exit(1); }
|
||||
#define s_stop(s, n) {exit(0);}
|
||||
#define z_abs(z) (cabs(Cd(z)))
|
||||
#define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
|
||||
#define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
|
||||
#define myexit_() break;
|
||||
#define mycycle_() continue;
|
||||
#define myceiling_(w) {ceil(w)}
|
||||
#define myhuge_(w) {HUGE_VAL}
|
||||
#define mymaxloc_(w,s,e,n) dmaxloc_(w,*(s),*(e),n)
|
||||
|
||||
/* procedure parameter types for -A and -C++ */
|
||||
|
||||
#define F2C_proc_par_types 1
|
||||
#ifdef __cplusplus
|
||||
typedef logical (*L_fp)(...);
|
||||
#else
|
||||
typedef logical (*L_fp)();
|
||||
#endif
|
||||
|
||||
/* -- translated by f2c (version 20000121).
|
||||
You must link the resulting object file with the libraries:
|
||||
-lf2c -lm (in that order)
|
||||
*/
|
||||
|
||||
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static real c_b9 = 1.f;
|
||||
|
||||
/* > \brief \b SGEQRS */
|
||||
|
||||
/* =========== DOCUMENTATION =========== */
|
||||
|
||||
/* Online html documentation available at */
|
||||
/* http://www.netlib.org/lapack/explore-html/ */
|
||||
|
||||
/* Definition: */
|
||||
/* =========== */
|
||||
|
||||
/* SUBROUTINE SGEQRS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, */
|
||||
/* INFO ) */
|
||||
|
||||
/* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS */
|
||||
/* REAL A( LDA, * ), B( LDB, * ), TAU( * ), */
|
||||
/* $ WORK( LWORK ) */
|
||||
|
||||
|
||||
/* > \par Purpose: */
|
||||
/* ============= */
|
||||
/* > */
|
||||
/* > \verbatim */
|
||||
/* > */
|
||||
/* > Solve the least squares problem */
|
||||
/* > f2cmin || A*X - B || */
|
||||
/* > using the QR factorization */
|
||||
/* > A = Q*R */
|
||||
/* > computed by SGEQRF. */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Arguments: */
|
||||
/* ========== */
|
||||
|
||||
/* > \param[in] M */
|
||||
/* > \verbatim */
|
||||
/* > M is INTEGER */
|
||||
/* > The number of rows of the matrix A. M >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] N */
|
||||
/* > \verbatim */
|
||||
/* > N is INTEGER */
|
||||
/* > The number of columns of the matrix A. M >= N >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] NRHS */
|
||||
/* > \verbatim */
|
||||
/* > NRHS is INTEGER */
|
||||
/* > The number of columns of B. NRHS >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] A */
|
||||
/* > \verbatim */
|
||||
/* > A is REAL array, dimension (LDA,N) */
|
||||
/* > Details of the QR factorization of the original matrix A as */
|
||||
/* > returned by SGEQRF. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDA */
|
||||
/* > \verbatim */
|
||||
/* > LDA is INTEGER */
|
||||
/* > The leading dimension of the array A. LDA >= M. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] TAU */
|
||||
/* > \verbatim */
|
||||
/* > TAU is REAL array, dimension (N) */
|
||||
/* > Details of the orthogonal matrix Q. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in,out] B */
|
||||
/* > \verbatim */
|
||||
/* > B is REAL array, dimension (LDB,NRHS) */
|
||||
/* > On entry, the m-by-nrhs right hand side matrix B. */
|
||||
/* > On exit, the n-by-nrhs solution matrix X. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDB */
|
||||
/* > \verbatim */
|
||||
/* > LDB is INTEGER */
|
||||
/* > The leading dimension of the array B. LDB >= M. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] WORK */
|
||||
/* > \verbatim */
|
||||
/* > WORK is REAL array, dimension (LWORK) */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LWORK */
|
||||
/* > \verbatim */
|
||||
/* > LWORK is INTEGER */
|
||||
/* > The length of the array WORK. LWORK must be at least NRHS, */
|
||||
/* > and should be at least NRHS*NB, where NB is the block size */
|
||||
/* > for this environment. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] INFO */
|
||||
/* > \verbatim */
|
||||
/* > INFO is INTEGER */
|
||||
/* > = 0: successful exit */
|
||||
/* > < 0: if INFO = -i, the i-th argument had an illegal value */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Authors: */
|
||||
/* ======== */
|
||||
|
||||
/* > \author Univ. of Tennessee */
|
||||
/* > \author Univ. of California Berkeley */
|
||||
/* > \author Univ. of Colorado Denver */
|
||||
/* > \author NAG Ltd. */
|
||||
|
||||
/* > \ingroup single_lin */
|
||||
|
||||
/* ===================================================================== */
|
||||
/* Subroutine */ int sgeqrs_(integer *m, integer *n, integer *nrhs, real *a,
|
||||
integer *lda, real *tau, real *b, integer *ldb, real *work, integer *
|
||||
lwork, integer *info)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, b_dim1, b_offset, i__1;
|
||||
|
||||
/* Local variables */
|
||||
extern /* Subroutine */ int strsm_(char *, char *, char *, char *,
|
||||
integer *, integer *, real *, real *, integer *, real *, integer *
|
||||
), xerbla_(char *, integer *), sormqr_(char *, char *, integer *, integer *, integer *,
|
||||
real *, integer *, real *, real *, integer *, real *, integer *,
|
||||
integer *);
|
||||
|
||||
|
||||
/* -- LAPACK test routine -- */
|
||||
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
|
||||
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
|
||||
|
||||
|
||||
/* ===================================================================== */
|
||||
|
||||
|
||||
/* Test the input arguments. */
|
||||
|
||||
/* Parameter adjustments */
|
||||
a_dim1 = *lda;
|
||||
a_offset = 1 + a_dim1 * 1;
|
||||
a -= a_offset;
|
||||
--tau;
|
||||
b_dim1 = *ldb;
|
||||
b_offset = 1 + b_dim1 * 1;
|
||||
b -= b_offset;
|
||||
--work;
|
||||
|
||||
/* Function Body */
|
||||
*info = 0;
|
||||
if (*m < 0) {
|
||||
*info = -1;
|
||||
} else if (*n < 0 || *n > *m) {
|
||||
*info = -2;
|
||||
} else if (*nrhs < 0) {
|
||||
*info = -3;
|
||||
} else if (*lda < f2cmax(1,*m)) {
|
||||
*info = -5;
|
||||
} else if (*ldb < f2cmax(1,*m)) {
|
||||
*info = -8;
|
||||
} else if (*lwork < 1 || *lwork < *nrhs && *m > 0 && *n > 0) {
|
||||
*info = -10;
|
||||
}
|
||||
if (*info != 0) {
|
||||
i__1 = -(*info);
|
||||
xerbla_("SGEQRS", &i__1);
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Quick return if possible */
|
||||
|
||||
if (*n == 0 || *nrhs == 0 || *m == 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* B := Q' * B */
|
||||
|
||||
sormqr_("Left", "Transpose", m, nrhs, n, &a[a_offset], lda, &tau[1], &b[
|
||||
b_offset], ldb, &work[1], lwork, info);
|
||||
|
||||
/* Solve R*X = B(1:n,:) */
|
||||
|
||||
strsm_("Left", "Upper", "No transpose", "Non-unit", n, nrhs, &c_b9, &a[
|
||||
a_offset], lda, &b[b_offset], ldb);
|
||||
|
||||
return 0;
|
||||
|
||||
/* End of SGEQRS */
|
||||
|
||||
} /* sgeqrs_ */
|
||||
|
|
@ -0,0 +1,481 @@
|
|||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <stdio.h>
|
||||
#include <complex.h>
|
||||
#ifdef complex
|
||||
#undef complex
|
||||
#endif
|
||||
#ifdef I
|
||||
#undef I
|
||||
#endif
|
||||
|
||||
#if defined(_WIN64)
|
||||
typedef long long BLASLONG;
|
||||
typedef unsigned long long BLASULONG;
|
||||
#else
|
||||
typedef long BLASLONG;
|
||||
typedef unsigned long BLASULONG;
|
||||
#endif
|
||||
|
||||
#ifdef LAPACK_ILP64
|
||||
typedef BLASLONG blasint;
|
||||
#if defined(_WIN64)
|
||||
#define blasabs(x) llabs(x)
|
||||
#else
|
||||
#define blasabs(x) labs(x)
|
||||
#endif
|
||||
#else
|
||||
typedef int blasint;
|
||||
#define blasabs(x) abs(x)
|
||||
#endif
|
||||
|
||||
typedef blasint integer;
|
||||
|
||||
typedef unsigned int uinteger;
|
||||
typedef char *address;
|
||||
typedef short int shortint;
|
||||
typedef float real;
|
||||
typedef double doublereal;
|
||||
typedef struct { real r, i; } complex;
|
||||
typedef struct { doublereal r, i; } doublecomplex;
|
||||
#ifdef _MSC_VER
|
||||
static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
|
||||
static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
|
||||
static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
|
||||
static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
|
||||
#else
|
||||
static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
|
||||
static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
|
||||
#endif
|
||||
#define pCf(z) (*_pCf(z))
|
||||
#define pCd(z) (*_pCd(z))
|
||||
typedef int logical;
|
||||
typedef short int shortlogical;
|
||||
typedef char logical1;
|
||||
typedef char integer1;
|
||||
|
||||
#define TRUE_ (1)
|
||||
#define FALSE_ (0)
|
||||
|
||||
/* Extern is for use with -E */
|
||||
#ifndef Extern
|
||||
#define Extern extern
|
||||
#endif
|
||||
|
||||
/* I/O stuff */
|
||||
|
||||
typedef int flag;
|
||||
typedef int ftnlen;
|
||||
typedef int ftnint;
|
||||
|
||||
/*external read, write*/
|
||||
typedef struct
|
||||
{ flag cierr;
|
||||
ftnint ciunit;
|
||||
flag ciend;
|
||||
char *cifmt;
|
||||
ftnint cirec;
|
||||
} cilist;
|
||||
|
||||
/*internal read, write*/
|
||||
typedef struct
|
||||
{ flag icierr;
|
||||
char *iciunit;
|
||||
flag iciend;
|
||||
char *icifmt;
|
||||
ftnint icirlen;
|
||||
ftnint icirnum;
|
||||
} icilist;
|
||||
|
||||
/*open*/
|
||||
typedef struct
|
||||
{ flag oerr;
|
||||
ftnint ounit;
|
||||
char *ofnm;
|
||||
ftnlen ofnmlen;
|
||||
char *osta;
|
||||
char *oacc;
|
||||
char *ofm;
|
||||
ftnint orl;
|
||||
char *oblnk;
|
||||
} olist;
|
||||
|
||||
/*close*/
|
||||
typedef struct
|
||||
{ flag cerr;
|
||||
ftnint cunit;
|
||||
char *csta;
|
||||
} cllist;
|
||||
|
||||
/*rewind, backspace, endfile*/
|
||||
typedef struct
|
||||
{ flag aerr;
|
||||
ftnint aunit;
|
||||
} alist;
|
||||
|
||||
/* inquire */
|
||||
typedef struct
|
||||
{ flag inerr;
|
||||
ftnint inunit;
|
||||
char *infile;
|
||||
ftnlen infilen;
|
||||
ftnint *inex; /*parameters in standard's order*/
|
||||
ftnint *inopen;
|
||||
ftnint *innum;
|
||||
ftnint *innamed;
|
||||
char *inname;
|
||||
ftnlen innamlen;
|
||||
char *inacc;
|
||||
ftnlen inacclen;
|
||||
char *inseq;
|
||||
ftnlen inseqlen;
|
||||
char *indir;
|
||||
ftnlen indirlen;
|
||||
char *infmt;
|
||||
ftnlen infmtlen;
|
||||
char *inform;
|
||||
ftnint informlen;
|
||||
char *inunf;
|
||||
ftnlen inunflen;
|
||||
ftnint *inrecl;
|
||||
ftnint *innrec;
|
||||
char *inblank;
|
||||
ftnlen inblanklen;
|
||||
} inlist;
|
||||
|
||||
#define VOID void
|
||||
|
||||
union Multitype { /* for multiple entry points */
|
||||
integer1 g;
|
||||
shortint h;
|
||||
integer i;
|
||||
/* longint j; */
|
||||
real r;
|
||||
doublereal d;
|
||||
complex c;
|
||||
doublecomplex z;
|
||||
};
|
||||
|
||||
typedef union Multitype Multitype;
|
||||
|
||||
struct Vardesc { /* for Namelist */
|
||||
char *name;
|
||||
char *addr;
|
||||
ftnlen *dims;
|
||||
int type;
|
||||
};
|
||||
typedef struct Vardesc Vardesc;
|
||||
|
||||
struct Namelist {
|
||||
char *name;
|
||||
Vardesc **vars;
|
||||
int nvars;
|
||||
};
|
||||
typedef struct Namelist Namelist;
|
||||
|
||||
#define abs(x) ((x) >= 0 ? (x) : -(x))
|
||||
#define dabs(x) (fabs(x))
|
||||
#define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
|
||||
#define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
|
||||
#define dmin(a,b) (f2cmin(a,b))
|
||||
#define dmax(a,b) (f2cmax(a,b))
|
||||
#define bit_test(a,b) ((a) >> (b) & 1)
|
||||
#define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
|
||||
#define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
|
||||
|
||||
#define abort_() { sig_die("Fortran abort routine called", 1); }
|
||||
#define c_abs(z) (cabsf(Cf(z)))
|
||||
#define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
|
||||
#ifdef _MSC_VER
|
||||
#define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
|
||||
#define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/Cd(b)._Val[1]);}
|
||||
#else
|
||||
#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
|
||||
#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
|
||||
#endif
|
||||
#define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
|
||||
#define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
|
||||
#define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
|
||||
//#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
|
||||
#define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
|
||||
#define d_abs(x) (fabs(*(x)))
|
||||
#define d_acos(x) (acos(*(x)))
|
||||
#define d_asin(x) (asin(*(x)))
|
||||
#define d_atan(x) (atan(*(x)))
|
||||
#define d_atn2(x, y) (atan2(*(x),*(y)))
|
||||
#define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
|
||||
#define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
|
||||
#define d_cos(x) (cos(*(x)))
|
||||
#define d_cosh(x) (cosh(*(x)))
|
||||
#define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
|
||||
#define d_exp(x) (exp(*(x)))
|
||||
#define d_imag(z) (cimag(Cd(z)))
|
||||
#define r_imag(z) (cimagf(Cf(z)))
|
||||
#define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define d_log(x) (log(*(x)))
|
||||
#define d_mod(x, y) (fmod(*(x), *(y)))
|
||||
#define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
|
||||
#define d_nint(x) u_nint(*(x))
|
||||
#define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
|
||||
#define d_sign(a,b) u_sign(*(a),*(b))
|
||||
#define r_sign(a,b) u_sign(*(a),*(b))
|
||||
#define d_sin(x) (sin(*(x)))
|
||||
#define d_sinh(x) (sinh(*(x)))
|
||||
#define d_sqrt(x) (sqrt(*(x)))
|
||||
#define d_tan(x) (tan(*(x)))
|
||||
#define d_tanh(x) (tanh(*(x)))
|
||||
#define i_abs(x) abs(*(x))
|
||||
#define i_dnnt(x) ((integer)u_nint(*(x)))
|
||||
#define i_len(s, n) (n)
|
||||
#define i_nint(x) ((integer)u_nint(*(x)))
|
||||
#define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
|
||||
#define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
|
||||
#define pow_si(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_ri(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_di(B,E) dpow_ui(*(B),*(E))
|
||||
#define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
|
||||
#define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
|
||||
#define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
|
||||
#define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
|
||||
#define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
|
||||
#define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
|
||||
#define sig_die(s, kill) { exit(1); }
|
||||
#define s_stop(s, n) {exit(0);}
|
||||
#define z_abs(z) (cabs(Cd(z)))
|
||||
#define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
|
||||
#define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
|
||||
#define myexit_() break;
|
||||
#define mycycle_() continue;
|
||||
#define myceiling_(w) {ceil(w)}
|
||||
#define myhuge_(w) {HUGE_VAL}
|
||||
#define mymaxloc_(w,s,e,n) dmaxloc_(w,*(s),*(e),n)
|
||||
|
||||
/* procedure parameter types for -A and -C++ */
|
||||
|
||||
#define F2C_proc_par_types 1
|
||||
#ifdef __cplusplus
|
||||
typedef logical (*L_fp)(...);
|
||||
#else
|
||||
typedef logical (*L_fp)();
|
||||
#endif
|
||||
|
||||
/* -- translated by f2c (version 20000121).
|
||||
You must link the resulting object file with the libraries:
|
||||
-lf2c -lm (in that order)
|
||||
*/
|
||||
|
||||
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static doublecomplex c_b1 = {0.,0.};
|
||||
static doublecomplex c_b2 = {1.,0.};
|
||||
|
||||
/* > \brief \b ZGELQS */
|
||||
|
||||
/* =========== DOCUMENTATION =========== */
|
||||
|
||||
/* Online html documentation available at */
|
||||
/* http://www.netlib.org/lapack/explore-html/ */
|
||||
|
||||
/* Definition: */
|
||||
/* =========== */
|
||||
|
||||
/* SUBROUTINE ZGELQS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, */
|
||||
/* INFO ) */
|
||||
|
||||
/* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS */
|
||||
/* COMPLEX*16 A( LDA, * ), B( LDB, * ), TAU( * ), */
|
||||
/* $ WORK( LWORK ) */
|
||||
|
||||
|
||||
/* > \par Purpose: */
|
||||
/* ============= */
|
||||
/* > */
|
||||
/* > \verbatim */
|
||||
/* > */
|
||||
/* > Compute a minimum-norm solution */
|
||||
/* > f2cmin || A*X - B || */
|
||||
/* > using the LQ factorization */
|
||||
/* > A = L*Q */
|
||||
/* > computed by ZGELQF. */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Arguments: */
|
||||
/* ========== */
|
||||
|
||||
/* > \param[in] M */
|
||||
/* > \verbatim */
|
||||
/* > M is INTEGER */
|
||||
/* > The number of rows of the matrix A. M >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] N */
|
||||
/* > \verbatim */
|
||||
/* > N is INTEGER */
|
||||
/* > The number of columns of the matrix A. N >= M >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] NRHS */
|
||||
/* > \verbatim */
|
||||
/* > NRHS is INTEGER */
|
||||
/* > The number of columns of B. NRHS >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] A */
|
||||
/* > \verbatim */
|
||||
/* > A is COMPLEX*16 array, dimension (LDA,N) */
|
||||
/* > Details of the LQ factorization of the original matrix A as */
|
||||
/* > returned by ZGELQF. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDA */
|
||||
/* > \verbatim */
|
||||
/* > LDA is INTEGER */
|
||||
/* > The leading dimension of the array A. LDA >= M. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] TAU */
|
||||
/* > \verbatim */
|
||||
/* > TAU is COMPLEX*16 array, dimension (M) */
|
||||
/* > Details of the orthogonal matrix Q. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in,out] B */
|
||||
/* > \verbatim */
|
||||
/* > B is COMPLEX*16 array, dimension (LDB,NRHS) */
|
||||
/* > On entry, the m-by-nrhs right hand side matrix B. */
|
||||
/* > On exit, the n-by-nrhs solution matrix X. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDB */
|
||||
/* > \verbatim */
|
||||
/* > LDB is INTEGER */
|
||||
/* > The leading dimension of the array B. LDB >= N. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] WORK */
|
||||
/* > \verbatim */
|
||||
/* > WORK is COMPLEX*16 array, dimension (LWORK) */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LWORK */
|
||||
/* > \verbatim */
|
||||
/* > LWORK is INTEGER */
|
||||
/* > The length of the array WORK. LWORK must be at least NRHS, */
|
||||
/* > and should be at least NRHS*NB, where NB is the block size */
|
||||
/* > for this environment. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] INFO */
|
||||
/* > \verbatim */
|
||||
/* > INFO is INTEGER */
|
||||
/* > = 0: successful exit */
|
||||
/* > < 0: if INFO = -i, the i-th argument had an illegal value */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Authors: */
|
||||
/* ======== */
|
||||
|
||||
/* > \author Univ. of Tennessee */
|
||||
/* > \author Univ. of California Berkeley */
|
||||
/* > \author Univ. of Colorado Denver */
|
||||
/* > \author NAG Ltd. */
|
||||
|
||||
/* > \ingroup complex16_lin */
|
||||
|
||||
/* ===================================================================== */
|
||||
/* Subroutine */ int zgelqs_(integer *m, integer *n, integer *nrhs,
|
||||
doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *b,
|
||||
integer *ldb, doublecomplex *work, integer *lwork, integer *info)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, b_dim1, b_offset, i__1;
|
||||
|
||||
/* Local variables */
|
||||
extern /* Subroutine */ int ztrsm_(char *, char *, char *, char *,
|
||||
integer *, integer *, doublecomplex *, doublecomplex *, integer *,
|
||||
doublecomplex *, integer *),
|
||||
xerbla_(char *, integer *), zlaset_(char *, integer *,
|
||||
integer *, doublecomplex *, doublecomplex *, doublecomplex *,
|
||||
integer *), zunmlq_(char *, char *, integer *, integer *,
|
||||
integer *, doublecomplex *, integer *, doublecomplex *,
|
||||
doublecomplex *, integer *, doublecomplex *, integer *, integer *);
|
||||
|
||||
|
||||
/* -- LAPACK test routine -- */
|
||||
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
|
||||
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
|
||||
|
||||
|
||||
/* ===================================================================== */
|
||||
|
||||
|
||||
/* Test the input parameters. */
|
||||
|
||||
/* Parameter adjustments */
|
||||
a_dim1 = *lda;
|
||||
a_offset = 1 + a_dim1 * 1;
|
||||
a -= a_offset;
|
||||
--tau;
|
||||
b_dim1 = *ldb;
|
||||
b_offset = 1 + b_dim1 * 1;
|
||||
b -= b_offset;
|
||||
--work;
|
||||
|
||||
/* Function Body */
|
||||
*info = 0;
|
||||
if (*m < 0) {
|
||||
*info = -1;
|
||||
} else if (*n < 0 || *m > *n) {
|
||||
*info = -2;
|
||||
} else if (*nrhs < 0) {
|
||||
*info = -3;
|
||||
} else if (*lda < f2cmax(1,*m)) {
|
||||
*info = -5;
|
||||
} else if (*ldb < f2cmax(1,*n)) {
|
||||
*info = -8;
|
||||
} else if (*lwork < 1 || *lwork < *nrhs && *m > 0 && *n > 0) {
|
||||
*info = -10;
|
||||
}
|
||||
if (*info != 0) {
|
||||
i__1 = -(*info);
|
||||
xerbla_("ZGELQS", &i__1);
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Quick return if possible */
|
||||
|
||||
if (*n == 0 || *nrhs == 0 || *m == 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Solve L*X = B(1:m,:) */
|
||||
|
||||
ztrsm_("Left", "Lower", "No transpose", "Non-unit", m, nrhs, &c_b2, &a[
|
||||
a_offset], lda, &b[b_offset], ldb);
|
||||
|
||||
/* Set B(m+1:n,:) to zero */
|
||||
|
||||
if (*m < *n) {
|
||||
i__1 = *n - *m;
|
||||
zlaset_("Full", &i__1, nrhs, &c_b1, &c_b1, &b[*m + 1 + b_dim1], ldb);
|
||||
}
|
||||
|
||||
/* B := Q' * B */
|
||||
|
||||
zunmlq_("Left", "Conjugate transpose", n, nrhs, m, &a[a_offset], lda, &
|
||||
tau[1], &b[b_offset], ldb, &work[1], lwork, info);
|
||||
|
||||
return 0;
|
||||
|
||||
/* End of ZGELQS */
|
||||
|
||||
} /* zgelqs_ */
|
||||
|
|
@ -0,0 +1,472 @@
|
|||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <stdio.h>
|
||||
#include <complex.h>
|
||||
#ifdef complex
|
||||
#undef complex
|
||||
#endif
|
||||
#ifdef I
|
||||
#undef I
|
||||
#endif
|
||||
|
||||
#if defined(_WIN64)
|
||||
typedef long long BLASLONG;
|
||||
typedef unsigned long long BLASULONG;
|
||||
#else
|
||||
typedef long BLASLONG;
|
||||
typedef unsigned long BLASULONG;
|
||||
#endif
|
||||
|
||||
#ifdef LAPACK_ILP64
|
||||
typedef BLASLONG blasint;
|
||||
#if defined(_WIN64)
|
||||
#define blasabs(x) llabs(x)
|
||||
#else
|
||||
#define blasabs(x) labs(x)
|
||||
#endif
|
||||
#else
|
||||
typedef int blasint;
|
||||
#define blasabs(x) abs(x)
|
||||
#endif
|
||||
|
||||
typedef blasint integer;
|
||||
|
||||
typedef unsigned int uinteger;
|
||||
typedef char *address;
|
||||
typedef short int shortint;
|
||||
typedef float real;
|
||||
typedef double doublereal;
|
||||
typedef struct { real r, i; } complex;
|
||||
typedef struct { doublereal r, i; } doublecomplex;
|
||||
#ifdef _MSC_VER
|
||||
static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
|
||||
static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
|
||||
static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
|
||||
static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
|
||||
#else
|
||||
static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
|
||||
static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
|
||||
static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
|
||||
#endif
|
||||
#define pCf(z) (*_pCf(z))
|
||||
#define pCd(z) (*_pCd(z))
|
||||
typedef int logical;
|
||||
typedef short int shortlogical;
|
||||
typedef char logical1;
|
||||
typedef char integer1;
|
||||
|
||||
#define TRUE_ (1)
|
||||
#define FALSE_ (0)
|
||||
|
||||
/* Extern is for use with -E */
|
||||
#ifndef Extern
|
||||
#define Extern extern
|
||||
#endif
|
||||
|
||||
/* I/O stuff */
|
||||
|
||||
typedef int flag;
|
||||
typedef int ftnlen;
|
||||
typedef int ftnint;
|
||||
|
||||
/*external read, write*/
|
||||
typedef struct
|
||||
{ flag cierr;
|
||||
ftnint ciunit;
|
||||
flag ciend;
|
||||
char *cifmt;
|
||||
ftnint cirec;
|
||||
} cilist;
|
||||
|
||||
/*internal read, write*/
|
||||
typedef struct
|
||||
{ flag icierr;
|
||||
char *iciunit;
|
||||
flag iciend;
|
||||
char *icifmt;
|
||||
ftnint icirlen;
|
||||
ftnint icirnum;
|
||||
} icilist;
|
||||
|
||||
/*open*/
|
||||
typedef struct
|
||||
{ flag oerr;
|
||||
ftnint ounit;
|
||||
char *ofnm;
|
||||
ftnlen ofnmlen;
|
||||
char *osta;
|
||||
char *oacc;
|
||||
char *ofm;
|
||||
ftnint orl;
|
||||
char *oblnk;
|
||||
} olist;
|
||||
|
||||
/*close*/
|
||||
typedef struct
|
||||
{ flag cerr;
|
||||
ftnint cunit;
|
||||
char *csta;
|
||||
} cllist;
|
||||
|
||||
/*rewind, backspace, endfile*/
|
||||
typedef struct
|
||||
{ flag aerr;
|
||||
ftnint aunit;
|
||||
} alist;
|
||||
|
||||
/* inquire */
|
||||
typedef struct
|
||||
{ flag inerr;
|
||||
ftnint inunit;
|
||||
char *infile;
|
||||
ftnlen infilen;
|
||||
ftnint *inex; /*parameters in standard's order*/
|
||||
ftnint *inopen;
|
||||
ftnint *innum;
|
||||
ftnint *innamed;
|
||||
char *inname;
|
||||
ftnlen innamlen;
|
||||
char *inacc;
|
||||
ftnlen inacclen;
|
||||
char *inseq;
|
||||
ftnlen inseqlen;
|
||||
char *indir;
|
||||
ftnlen indirlen;
|
||||
char *infmt;
|
||||
ftnlen infmtlen;
|
||||
char *inform;
|
||||
ftnint informlen;
|
||||
char *inunf;
|
||||
ftnlen inunflen;
|
||||
ftnint *inrecl;
|
||||
ftnint *innrec;
|
||||
char *inblank;
|
||||
ftnlen inblanklen;
|
||||
} inlist;
|
||||
|
||||
#define VOID void
|
||||
|
||||
union Multitype { /* for multiple entry points */
|
||||
integer1 g;
|
||||
shortint h;
|
||||
integer i;
|
||||
/* longint j; */
|
||||
real r;
|
||||
doublereal d;
|
||||
complex c;
|
||||
doublecomplex z;
|
||||
};
|
||||
|
||||
typedef union Multitype Multitype;
|
||||
|
||||
struct Vardesc { /* for Namelist */
|
||||
char *name;
|
||||
char *addr;
|
||||
ftnlen *dims;
|
||||
int type;
|
||||
};
|
||||
typedef struct Vardesc Vardesc;
|
||||
|
||||
struct Namelist {
|
||||
char *name;
|
||||
Vardesc **vars;
|
||||
int nvars;
|
||||
};
|
||||
typedef struct Namelist Namelist;
|
||||
|
||||
#define abs(x) ((x) >= 0 ? (x) : -(x))
|
||||
#define dabs(x) (fabs(x))
|
||||
#define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
|
||||
#define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
|
||||
#define dmin(a,b) (f2cmin(a,b))
|
||||
#define dmax(a,b) (f2cmax(a,b))
|
||||
#define bit_test(a,b) ((a) >> (b) & 1)
|
||||
#define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
|
||||
#define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
|
||||
|
||||
#define abort_() { sig_die("Fortran abort routine called", 1); }
|
||||
#define c_abs(z) (cabsf(Cf(z)))
|
||||
#define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
|
||||
#ifdef _MSC_VER
|
||||
#define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
|
||||
#define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/Cd(b)._Val[1]);}
|
||||
#else
|
||||
#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
|
||||
#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
|
||||
#endif
|
||||
#define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
|
||||
#define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
|
||||
#define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
|
||||
//#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
|
||||
#define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
|
||||
#define d_abs(x) (fabs(*(x)))
|
||||
#define d_acos(x) (acos(*(x)))
|
||||
#define d_asin(x) (asin(*(x)))
|
||||
#define d_atan(x) (atan(*(x)))
|
||||
#define d_atn2(x, y) (atan2(*(x),*(y)))
|
||||
#define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
|
||||
#define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
|
||||
#define d_cos(x) (cos(*(x)))
|
||||
#define d_cosh(x) (cosh(*(x)))
|
||||
#define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
|
||||
#define d_exp(x) (exp(*(x)))
|
||||
#define d_imag(z) (cimag(Cd(z)))
|
||||
#define r_imag(z) (cimagf(Cf(z)))
|
||||
#define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
|
||||
#define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
|
||||
#define d_log(x) (log(*(x)))
|
||||
#define d_mod(x, y) (fmod(*(x), *(y)))
|
||||
#define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
|
||||
#define d_nint(x) u_nint(*(x))
|
||||
#define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
|
||||
#define d_sign(a,b) u_sign(*(a),*(b))
|
||||
#define r_sign(a,b) u_sign(*(a),*(b))
|
||||
#define d_sin(x) (sin(*(x)))
|
||||
#define d_sinh(x) (sinh(*(x)))
|
||||
#define d_sqrt(x) (sqrt(*(x)))
|
||||
#define d_tan(x) (tan(*(x)))
|
||||
#define d_tanh(x) (tanh(*(x)))
|
||||
#define i_abs(x) abs(*(x))
|
||||
#define i_dnnt(x) ((integer)u_nint(*(x)))
|
||||
#define i_len(s, n) (n)
|
||||
#define i_nint(x) ((integer)u_nint(*(x)))
|
||||
#define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
|
||||
#define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
|
||||
#define pow_si(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_ri(B,E) spow_ui(*(B),*(E))
|
||||
#define pow_di(B,E) dpow_ui(*(B),*(E))
|
||||
#define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
|
||||
#define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
|
||||
#define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
|
||||
#define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
|
||||
#define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
|
||||
#define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
|
||||
#define sig_die(s, kill) { exit(1); }
|
||||
#define s_stop(s, n) {exit(0);}
|
||||
#define z_abs(z) (cabs(Cd(z)))
|
||||
#define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
|
||||
#define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
|
||||
#define myexit_() break;
|
||||
#define mycycle_() continue;
|
||||
#define myceiling_(w) {ceil(w)}
|
||||
#define myhuge_(w) {HUGE_VAL}
|
||||
#define mymaxloc_(w,s,e,n) dmaxloc_(w,*(s),*(e),n)
|
||||
|
||||
/* procedure parameter types for -A and -C++ */
|
||||
|
||||
#define F2C_proc_par_types 1
|
||||
#ifdef __cplusplus
|
||||
typedef logical (*L_fp)(...);
|
||||
#else
|
||||
typedef logical (*L_fp)();
|
||||
#endif
|
||||
|
||||
/* -- translated by f2c (version 20000121).
|
||||
You must link the resulting object file with the libraries:
|
||||
-lf2c -lm (in that order)
|
||||
*/
|
||||
|
||||
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static doublecomplex c_b1 = {1.,0.};
|
||||
|
||||
/* > \brief \b ZGEQRS */
|
||||
|
||||
/* =========== DOCUMENTATION =========== */
|
||||
|
||||
/* Online html documentation available at */
|
||||
/* http://www.netlib.org/lapack/explore-html/ */
|
||||
|
||||
/* Definition: */
|
||||
/* =========== */
|
||||
|
||||
/* SUBROUTINE ZGEQRS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, */
|
||||
/* INFO ) */
|
||||
|
||||
/* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS */
|
||||
/* COMPLEX*16 A( LDA, * ), B( LDB, * ), TAU( * ), */
|
||||
/* $ WORK( LWORK ) */
|
||||
|
||||
|
||||
/* > \par Purpose: */
|
||||
/* ============= */
|
||||
/* > */
|
||||
/* > \verbatim */
|
||||
/* > */
|
||||
/* > Solve the least squares problem */
|
||||
/* > f2cmin || A*X - B || */
|
||||
/* > using the QR factorization */
|
||||
/* > A = Q*R */
|
||||
/* > computed by ZGEQRF. */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Arguments: */
|
||||
/* ========== */
|
||||
|
||||
/* > \param[in] M */
|
||||
/* > \verbatim */
|
||||
/* > M is INTEGER */
|
||||
/* > The number of rows of the matrix A. M >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] N */
|
||||
/* > \verbatim */
|
||||
/* > N is INTEGER */
|
||||
/* > The number of columns of the matrix A. M >= N >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] NRHS */
|
||||
/* > \verbatim */
|
||||
/* > NRHS is INTEGER */
|
||||
/* > The number of columns of B. NRHS >= 0. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] A */
|
||||
/* > \verbatim */
|
||||
/* > A is COMPLEX*16 array, dimension (LDA,N) */
|
||||
/* > Details of the QR factorization of the original matrix A as */
|
||||
/* > returned by ZGEQRF. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDA */
|
||||
/* > \verbatim */
|
||||
/* > LDA is INTEGER */
|
||||
/* > The leading dimension of the array A. LDA >= M. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] TAU */
|
||||
/* > \verbatim */
|
||||
/* > TAU is COMPLEX*16 array, dimension (N) */
|
||||
/* > Details of the orthogonal matrix Q. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in,out] B */
|
||||
/* > \verbatim */
|
||||
/* > B is COMPLEX*16 array, dimension (LDB,NRHS) */
|
||||
/* > On entry, the m-by-nrhs right hand side matrix B. */
|
||||
/* > On exit, the n-by-nrhs solution matrix X. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LDB */
|
||||
/* > \verbatim */
|
||||
/* > LDB is INTEGER */
|
||||
/* > The leading dimension of the array B. LDB >= M. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] WORK */
|
||||
/* > \verbatim */
|
||||
/* > WORK is COMPLEX*16 array, dimension (LWORK) */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[in] LWORK */
|
||||
/* > \verbatim */
|
||||
/* > LWORK is INTEGER */
|
||||
/* > The length of the array WORK. LWORK must be at least NRHS, */
|
||||
/* > and should be at least NRHS*NB, where NB is the block size */
|
||||
/* > for this environment. */
|
||||
/* > \endverbatim */
|
||||
/* > */
|
||||
/* > \param[out] INFO */
|
||||
/* > \verbatim */
|
||||
/* > INFO is INTEGER */
|
||||
/* > = 0: successful exit */
|
||||
/* > < 0: if INFO = -i, the i-th argument had an illegal value */
|
||||
/* > \endverbatim */
|
||||
|
||||
/* Authors: */
|
||||
/* ======== */
|
||||
|
||||
/* > \author Univ. of Tennessee */
|
||||
/* > \author Univ. of California Berkeley */
|
||||
/* > \author Univ. of Colorado Denver */
|
||||
/* > \author NAG Ltd. */
|
||||
|
||||
/* > \ingroup complex16_lin */
|
||||
|
||||
/* ===================================================================== */
|
||||
/* Subroutine */ int zgeqrs_(integer *m, integer *n, integer *nrhs,
|
||||
doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *b,
|
||||
integer *ldb, doublecomplex *work, integer *lwork, integer *info)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, b_dim1, b_offset, i__1;
|
||||
|
||||
/* Local variables */
|
||||
extern /* Subroutine */ int ztrsm_(char *, char *, char *, char *,
|
||||
integer *, integer *, doublecomplex *, doublecomplex *, integer *,
|
||||
doublecomplex *, integer *),
|
||||
xerbla_(char *, integer *), zunmqr_(char *, char *,
|
||||
integer *, integer *, integer *, doublecomplex *, integer *,
|
||||
doublecomplex *, doublecomplex *, integer *, doublecomplex *,
|
||||
integer *, integer *);
|
||||
|
||||
|
||||
/* -- LAPACK test routine -- */
|
||||
/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
|
||||
/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
|
||||
|
||||
|
||||
/* ===================================================================== */
|
||||
|
||||
|
||||
/* Test the input arguments. */
|
||||
|
||||
/* Parameter adjustments */
|
||||
a_dim1 = *lda;
|
||||
a_offset = 1 + a_dim1 * 1;
|
||||
a -= a_offset;
|
||||
--tau;
|
||||
b_dim1 = *ldb;
|
||||
b_offset = 1 + b_dim1 * 1;
|
||||
b -= b_offset;
|
||||
--work;
|
||||
|
||||
/* Function Body */
|
||||
*info = 0;
|
||||
if (*m < 0) {
|
||||
*info = -1;
|
||||
} else if (*n < 0 || *n > *m) {
|
||||
*info = -2;
|
||||
} else if (*nrhs < 0) {
|
||||
*info = -3;
|
||||
} else if (*lda < f2cmax(1,*m)) {
|
||||
*info = -5;
|
||||
} else if (*ldb < f2cmax(1,*m)) {
|
||||
*info = -8;
|
||||
} else if (*lwork < 1 || *lwork < *nrhs && *m > 0 && *n > 0) {
|
||||
*info = -10;
|
||||
}
|
||||
if (*info != 0) {
|
||||
i__1 = -(*info);
|
||||
xerbla_("ZGEQRS", &i__1);
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Quick return if possible */
|
||||
|
||||
if (*n == 0 || *nrhs == 0 || *m == 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* B := Q' * B */
|
||||
|
||||
zunmqr_("Left", "Conjugate transpose", m, nrhs, n, &a[a_offset], lda, &
|
||||
tau[1], &b[b_offset], ldb, &work[1], lwork, info);
|
||||
|
||||
/* Solve R*X = B(1:n,:) */
|
||||
|
||||
ztrsm_("Left", "Upper", "No transpose", "Non-unit", n, nrhs, &c_b1, &a[
|
||||
a_offset], lda, &b[b_offset], ldb);
|
||||
|
||||
return 0;
|
||||
|
||||
/* End of ZGEQRS */
|
||||
|
||||
} /* zgeqrs_ */
|
||||
|
|
@ -544,26 +544,30 @@ endif
|
|||
ifeq ($(BUILD_COMPLEX),1)
|
||||
CDEPRECSRC = DEPRECATED/cgegs.o DEPRECATED/cgegv.o DEPRECATED/cgelsx.o \
|
||||
DEPRECATED/cgeqpf.o DEPRECATED/cggsvd.o DEPRECATED/cggsvp.o \
|
||||
DEPRECATED/clahrd.o DEPRECATED/clatzm.o DEPRECATED/ctzrqf.o
|
||||
DEPRECATED/clahrd.o DEPRECATED/clatzm.o DEPRECATED/ctzrqf.o \
|
||||
DEPRECATED/cgelqs.o DEPRECATED/cgeqrs.o
|
||||
endif
|
||||
|
||||
ifeq ($(BUILD_DOUBLE),1)
|
||||
DDEPRECSRC = \
|
||||
DEPRECATED/dgegs.o DEPRECATED/dgegv.o DEPRECATED/dgelsx.o \
|
||||
DEPRECATED/dgeqpf.o DEPRECATED/dggsvd.o DEPRECATED/dggsvp.o \
|
||||
DEPRECATED/dlahrd.o DEPRECATED/dlatzm.o DEPRECATED/dtzrqf.o
|
||||
DEPRECATED/dlahrd.o DEPRECATED/dlatzm.o DEPRECATED/dtzrqf.o \
|
||||
DEPRECATED/dgelqs.o DEPRECATED/dgeqrs.o
|
||||
endif
|
||||
ifeq ($(BUILD_SINGLE),1)
|
||||
SDEPRECSRC = \
|
||||
DEPRECATED/sgegs.o DEPRECATED/sgegv.o DEPRECATED/sgelsx.o \
|
||||
DEPRECATED/sgeqpf.o DEPRECATED/sggsvd.o DEPRECATED/sggsvp.o \
|
||||
DEPRECATED/slahrd.o DEPRECATED/slatzm.o DEPRECATED/stzrqf.o
|
||||
DEPRECATED/slahrd.o DEPRECATED/slatzm.o DEPRECATED/stzrqf.o \
|
||||
DEPRECATED/sgelqs.o DEPRECATED/sgeqrs.o
|
||||
endif
|
||||
ifeq ($(BUILD_COMPLEX16),1)
|
||||
ZDEPRECSRC = \
|
||||
DEPRECATED/zgegs.o DEPRECATED/zgegv.o DEPRECATED/zgelsx.o \
|
||||
DEPRECATED/zgeqpf.o DEPRECATED/zggsvd.o DEPRECATED/zggsvp.o \
|
||||
DEPRECATED/zlahrd.o DEPRECATED/zlatzm.o DEPRECATED/ztzrqf.o
|
||||
DEPRECATED/zlahrd.o DEPRECATED/zlatzm.o DEPRECATED/ztzrqf.o \
|
||||
DEPRECATED/zgelqs.o DEPRECATED/zgeqrs.o
|
||||
endif
|
||||
|
||||
# filter out optimized codes from OpenBLAS
|
||||
|
|
|
@ -20,7 +20,7 @@ set(SLINTST schkaa.F
|
|||
serrgt.f serrlq.f serrls.f
|
||||
serrps.f serrql.f serrqp.f serrqr.f
|
||||
serrrq.f serrtr.f serrtz.f
|
||||
sgbt01.f sgbt02.f sgbt05.f sgelqs.f sgeqls.f sgeqrs.f
|
||||
sgbt01.f sgbt02.f sgbt05.f sgeqls.f
|
||||
sgerqs.f sget01.f sget02.f
|
||||
sget03.f sget04.f sget06.f sget07.f sgtt01.f sgtt02.f
|
||||
sgtt05.f slaptm.f slarhs.f slatb4.f slatb5.f slattb.f slattp.f
|
||||
|
@ -70,7 +70,7 @@ set(CLINTST cchkaa.F
|
|||
cerrgt.f cerrlq.f
|
||||
cerrls.f cerrps.f cerrql.f cerrqp.f
|
||||
cerrqr.f cerrrq.f cerrtr.f cerrtz.f
|
||||
cgbt01.f cgbt02.f cgbt05.f cgelqs.f cgeqls.f cgeqrs.f
|
||||
cgbt01.f cgbt02.f cgbt05.f cgeqls.f
|
||||
cgerqs.f cget01.f cget02.f
|
||||
cget03.f cget04.f cget07.f cgtt01.f cgtt02.f
|
||||
cgtt05.f chet01.f chet01_rook.f chet01_3.f
|
||||
|
@ -121,7 +121,7 @@ set(DLINTST dchkaa.F
|
|||
derrgt.f derrlq.f derrls.f
|
||||
derrps.f derrql.f derrqp.f derrqr.f
|
||||
derrrq.f derrtr.f derrtz.f
|
||||
dgbt01.f dgbt02.f dgbt05.f dgelqs.f dgeqls.f dgeqrs.f
|
||||
dgbt01.f dgbt02.f dgbt05.f dgeqls.f
|
||||
dgerqs.f dget01.f dget02.f
|
||||
dget03.f dget04.f dget06.f dget07.f dgtt01.f dgtt02.f
|
||||
dgtt05.f dlaptm.f dlarhs.f dlatb4.f dlatb5.f dlattb.f dlattp.f
|
||||
|
@ -172,7 +172,7 @@ set(ZLINTST zchkaa.F
|
|||
zerrgt.f zerrlq.f
|
||||
zerrls.f zerrps.f zerrql.f zerrqp.f
|
||||
zerrqr.f zerrrq.f zerrtr.f zerrtz.f
|
||||
zgbt01.f zgbt02.f zgbt05.f zgelqs.f zgeqls.f zgeqrs.f
|
||||
zgbt01.f zgbt02.f zgbt05.f zgeqls.f
|
||||
zgerqs.f zget01.f zget02.f
|
||||
zget03.f zget04.f zget07.f zgtt01.f zgtt02.f
|
||||
zgtt05.f zhet01.f zhet01_rook.f zhet01_3.f
|
||||
|
|
|
@ -55,7 +55,7 @@ SLINTST = schkaa.o \
|
|||
serrgt.o serrlq.o serrls.o \
|
||||
serrps.o serrql.o serrqp.o serrqr.o \
|
||||
serrrq.o serrtr.o serrtz.o \
|
||||
sgbt01.o sgbt02.o sgbt05.o sgelqs.o sgeqls.o sgeqrs.o \
|
||||
sgbt01.o sgbt02.o sgbt05.o sgeqls.o \
|
||||
sgerqs.o sget01.o sget02.o \
|
||||
sget03.o sget04.o sget06.o sget07.o sgtt01.o sgtt02.o \
|
||||
sgtt05.o slaptm.o slarhs.o slatb4.o slatb5.o slattb.o slattp.o \
|
||||
|
@ -100,7 +100,7 @@ CLINTST = cchkaa.o \
|
|||
cerrgt.o cerrlq.o \
|
||||
cerrls.o cerrps.o cerrql.o cerrqp.o \
|
||||
cerrqr.o cerrrq.o cerrtr.o cerrtz.o \
|
||||
cgbt01.o cgbt02.o cgbt05.o cgelqs.o cgeqls.o cgeqrs.o \
|
||||
cgbt01.o cgbt02.o cgbt05.o cgeqls.o \
|
||||
cgerqs.o cget01.o cget02.o \
|
||||
cget03.o cget04.o cget07.o cgtt01.o cgtt02.o \
|
||||
cgtt05.o chet01.o chet01_rook.o chet01_3.o chet01_aa.o \
|
||||
|
@ -147,7 +147,7 @@ DLINTST = dchkaa.o \
|
|||
derrgt.o derrlq.o derrls.o \
|
||||
derrps.o derrql.o derrqp.o derrqr.o \
|
||||
derrrq.o derrtr.o derrtz.o \
|
||||
dgbt01.o dgbt02.o dgbt05.o dgelqs.o dgeqls.o dgeqrs.o \
|
||||
dgbt01.o dgbt02.o dgbt05.o dgeqls.o \
|
||||
dgerqs.o dget01.o dget02.o \
|
||||
dget03.o dget04.o dget06.o dget07.o dgtt01.o dgtt02.o \
|
||||
dgtt05.o dlaptm.o dlarhs.o dlatb4.o dlatb5.o dlattb.o dlattp.o \
|
||||
|
@ -192,7 +192,7 @@ ZLINTST = zchkaa.o \
|
|||
zerrgt.o zerrlq.o \
|
||||
zerrls.o zerrps.o zerrql.o zerrqp.o \
|
||||
zerrqr.o zerrrq.o zerrtr.o zerrtz.o \
|
||||
zgbt01.o zgbt02.o zgbt05.o zgelqs.o zgeqls.o zgeqrs.o \
|
||||
zgbt01.o zgbt02.o zgbt05.o zgeqls.o \
|
||||
zgerqs.o zget01.o zget02.o \
|
||||
zget03.o zget04.o zget07.o zgtt01.o zgtt02.o \
|
||||
zgtt05.o zhet01.o zhet01_rook.o zhet01_3.o zhet01_aa.o \
|
||||
|
|
|
@ -235,7 +235,7 @@
|
|||
REAL RESULT( NTESTS )
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, CERRLQ, CGELQS, CGET02,
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, CERRLQ, CGELS, CGET02,
|
||||
$ CLACPY, CLARHS, CLATB4, CLATMS, CLQT01, CLQT02,
|
||||
$ CLQT03, XLAENV
|
||||
* ..
|
||||
|
@ -370,7 +370,7 @@
|
|||
$ WORK, LWORK, RWORK, RESULT( 3 ) )
|
||||
NT = NT + 4
|
||||
*
|
||||
* If M>=N and K=N, call CGELQS to solve a system
|
||||
* If M<=N and K=M, call CGELS to solve a system
|
||||
* with NRHS right hand sides and compute the
|
||||
* residual.
|
||||
*
|
||||
|
@ -387,14 +387,20 @@
|
|||
*
|
||||
CALL CLACPY( 'Full', M, NRHS, B, LDA, X,
|
||||
$ LDA )
|
||||
SRNAMT = 'CGELQS'
|
||||
CALL CGELQS( M, N, NRHS, AF, LDA, TAU, X,
|
||||
$ LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from CGELQS.
|
||||
* Reset AF to the original matrix. CGELS
|
||||
* factors the matrix before solving the system.
|
||||
*
|
||||
CALL CLACPY( 'Full', M, N, A, LDA, AF, LDA )
|
||||
*
|
||||
SRNAMT = 'CGELS'
|
||||
CALL CGELS( 'No transpose', M, N, NRHS, AF,
|
||||
$ LDA, X, LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from CGELS.
|
||||
*
|
||||
IF( INFO.NE.0 )
|
||||
$ CALL ALAERH( PATH, 'CGELQS', INFO, 0, ' ',
|
||||
$ CALL ALAERH( PATH, 'CGELS', INFO, 0, 'N',
|
||||
$ M, N, NRHS, -1, NB, IMAT,
|
||||
$ NFAIL, NERRS, NOUT )
|
||||
*
|
||||
|
|
|
@ -244,7 +244,7 @@
|
|||
EXTERNAL CGENND
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, CERRQR, CGEQRS, CGET02,
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, CERRQR, CGELS, CGET02,
|
||||
$ CLACPY, CLARHS, CLATB4, CLATMS, CQRT01,
|
||||
$ CQRT01P, CQRT02, CQRT03, XLAENV
|
||||
* ..
|
||||
|
@ -371,7 +371,7 @@
|
|||
IF( .NOT. CGENND( M, N, AF, LDA ) )
|
||||
$ RESULT( 9 ) = 2*THRESH
|
||||
NT = NT + 1
|
||||
ELSE IF( M.GE.N ) THEN
|
||||
ELSE IF( M.GE.N ) THEN
|
||||
*
|
||||
* Test CUNGQR, using factorization
|
||||
* returned by CQRT01
|
||||
|
@ -388,7 +388,7 @@
|
|||
$ WORK, LWORK, RWORK, RESULT( 3 ) )
|
||||
NT = NT + 4
|
||||
*
|
||||
* If M>=N and K=N, call CGEQRS to solve a system
|
||||
* If M>=N and K=N, call CGELS to solve a system
|
||||
* with NRHS right hand sides and compute the
|
||||
* residual.
|
||||
*
|
||||
|
@ -405,14 +405,20 @@
|
|||
*
|
||||
CALL CLACPY( 'Full', M, NRHS, B, LDA, X,
|
||||
$ LDA )
|
||||
SRNAMT = 'CGEQRS'
|
||||
CALL CGEQRS( M, N, NRHS, AF, LDA, TAU, X,
|
||||
$ LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from CGEQRS.
|
||||
* Reset AF to the original matrix. CGELS
|
||||
* factors the matrix before solving the system.
|
||||
*
|
||||
CALL CLACPY( 'Full', M, N, A, LDA, AF, LDA )
|
||||
*
|
||||
SRNAMT = 'CGELS'
|
||||
CALL CGELS( 'No transpose', M, N, NRHS, AF,
|
||||
$ LDA, X, LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from CGELS.
|
||||
*
|
||||
IF( INFO.NE.0 )
|
||||
$ CALL ALAERH( PATH, 'CGEQRS', INFO, 0, ' ',
|
||||
$ CALL ALAERH( PATH, 'CGELS', INFO, 0, 'N',
|
||||
$ M, N, NRHS, -1, NB, IMAT,
|
||||
$ NFAIL, NERRS, NOUT )
|
||||
*
|
||||
|
|
|
@ -76,7 +76,7 @@
|
|||
$ W( NMAX ), X( NMAX )
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAESM, CGELQ2, CGELQF, CGELQS, CHKXER, CUNGL2,
|
||||
EXTERNAL ALAESM, CGELQ2, CGELQF, CHKXER, CUNGL2,
|
||||
$ CUNGLQ, CUNML2, CUNMLQ
|
||||
* ..
|
||||
* .. Scalars in Common ..
|
||||
|
@ -140,31 +140,6 @@
|
|||
CALL CGELQ2( 2, 1, A, 1, B, W, INFO )
|
||||
CALL CHKXER( 'CGELQ2', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* CGELQS
|
||||
*
|
||||
SRNAMT = 'CGELQS'
|
||||
INFOT = 1
|
||||
CALL CGELQS( -1, 0, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'CGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL CGELQS( 0, -1, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'CGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL CGELQS( 2, 1, 0, A, 2, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'CGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 3
|
||||
CALL CGELQS( 0, 0, -1, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'CGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 5
|
||||
CALL CGELQS( 2, 2, 0, A, 1, X, B, 2, W, 1, INFO )
|
||||
CALL CHKXER( 'CGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 8
|
||||
CALL CGELQS( 1, 2, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'CGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 10
|
||||
CALL CGELQS( 1, 1, 2, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'CGELQS', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* CUNGLQ
|
||||
*
|
||||
SRNAMT = 'CUNGLQ'
|
||||
|
|
|
@ -77,7 +77,7 @@
|
|||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAESM, CGEQR2, CGEQR2P, CGEQRF, CGEQRFP,
|
||||
$ CGEQRS, CHKXER, CUNG2R, CUNGQR, CUNM2R,
|
||||
$ CHKXER, CUNG2R, CUNGQR, CUNM2R,
|
||||
$ CUNMQR
|
||||
* ..
|
||||
* .. Scalars in Common ..
|
||||
|
@ -170,31 +170,6 @@
|
|||
CALL CGEQR2P( 2, 1, A, 1, B, W, INFO )
|
||||
CALL CHKXER( 'CGEQR2P', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* CGEQRS
|
||||
*
|
||||
SRNAMT = 'CGEQRS'
|
||||
INFOT = 1
|
||||
CALL CGEQRS( -1, 0, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'CGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL CGEQRS( 0, -1, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'CGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL CGEQRS( 1, 2, 0, A, 2, X, B, 2, W, 1, INFO )
|
||||
CALL CHKXER( 'CGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 3
|
||||
CALL CGEQRS( 0, 0, -1, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'CGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 5
|
||||
CALL CGEQRS( 2, 1, 0, A, 1, X, B, 2, W, 1, INFO )
|
||||
CALL CHKXER( 'CGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 8
|
||||
CALL CGEQRS( 2, 1, 0, A, 2, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'CGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 10
|
||||
CALL CGEQRS( 1, 1, 2, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'CGEQRS', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* CUNGQR
|
||||
*
|
||||
SRNAMT = 'CUNGQR'
|
||||
|
|
|
@ -235,7 +235,7 @@
|
|||
DOUBLE PRECISION RESULT( NTESTS )
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, DERRLQ, DGELQS, DGET02,
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, DERRLQ, DGELS, DGET02,
|
||||
$ DLACPY, DLARHS, DLATB4, DLATMS, DLQT01, DLQT02,
|
||||
$ DLQT03, XLAENV
|
||||
* ..
|
||||
|
@ -373,7 +373,7 @@
|
|||
$ WORK, LWORK, RWORK, RESULT( 3 ) )
|
||||
NT = NT + 4
|
||||
*
|
||||
* If M>=N and K=N, call DGELQS to solve a system
|
||||
* If M<=N and K=M, call DGELS to solve a system
|
||||
* with NRHS right hand sides and compute the
|
||||
* residual.
|
||||
*
|
||||
|
@ -390,14 +390,20 @@
|
|||
*
|
||||
CALL DLACPY( 'Full', M, NRHS, B, LDA, X,
|
||||
$ LDA )
|
||||
SRNAMT = 'DGELQS'
|
||||
CALL DGELQS( M, N, NRHS, AF, LDA, TAU, X,
|
||||
$ LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from DGELQS.
|
||||
* Reset AF to the original matrix. DGELS
|
||||
* factors the matrix before solving the system.
|
||||
*
|
||||
CALL DLACPY( 'Full', M, N, A, LDA, AF, LDA )
|
||||
*
|
||||
SRNAMT = 'DGELS'
|
||||
CALL DGELS( 'No transpose', M, N, NRHS, AF,
|
||||
$ LDA, X, LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from DGELS.
|
||||
*
|
||||
IF( INFO.NE.0 )
|
||||
$ CALL ALAERH( PATH, 'DGELQS', INFO, 0, ' ',
|
||||
$ CALL ALAERH( PATH, 'DGELS', INFO, 0, 'N',
|
||||
$ M, N, NRHS, -1, NB, IMAT,
|
||||
$ NFAIL, NERRS, NOUT )
|
||||
*
|
||||
|
|
|
@ -244,7 +244,7 @@
|
|||
EXTERNAL DGENND
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, DERRQR, DGEQRS, DGET02,
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, DERRQR, DGELS, DGET02,
|
||||
$ DLACPY, DLARHS, DLATB4, DLATMS, DQRT01,
|
||||
$ DQRT01P, DQRT02, DQRT03, XLAENV
|
||||
* ..
|
||||
|
@ -372,7 +372,7 @@
|
|||
IF( .NOT. DGENND( M, N, AF, LDA ) )
|
||||
$ RESULT( 9 ) = 2*THRESH
|
||||
NT = NT + 1
|
||||
ELSE IF( M.GE.N ) THEN
|
||||
ELSE IF( M.GE.N ) THEN
|
||||
*
|
||||
* Test DORGQR, using factorization
|
||||
* returned by DQRT01
|
||||
|
@ -389,7 +389,7 @@
|
|||
$ WORK, LWORK, RWORK, RESULT( 3 ) )
|
||||
NT = NT + 4
|
||||
*
|
||||
* If M>=N and K=N, call DGEQRS to solve a system
|
||||
* If M>=N and K=N, call DGELS to solve a system
|
||||
* with NRHS right hand sides and compute the
|
||||
* residual.
|
||||
*
|
||||
|
@ -406,14 +406,20 @@
|
|||
*
|
||||
CALL DLACPY( 'Full', M, NRHS, B, LDA, X,
|
||||
$ LDA )
|
||||
SRNAMT = 'DGEQRS'
|
||||
CALL DGEQRS( M, N, NRHS, AF, LDA, TAU, X,
|
||||
$ LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from DGEQRS.
|
||||
* Reset AF. DGELS overwrites the matrix with
|
||||
* its factorization.
|
||||
*
|
||||
CALL DLACPY( 'Full', M, N, A, LDA, AF, LDA )
|
||||
*
|
||||
SRNAMT = 'DGELS'
|
||||
CALL DGELS( 'No transpose', M, N, NRHS, AF,
|
||||
$ LDA, X, LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from DGELS.
|
||||
*
|
||||
IF( INFO.NE.0 )
|
||||
$ CALL ALAERH( PATH, 'DGEQRS', INFO, 0, ' ',
|
||||
$ CALL ALAERH( PATH, 'DGELS', INFO, 0, 'N',
|
||||
$ M, N, NRHS, -1, NB, IMAT,
|
||||
$ NFAIL, NERRS, NOUT )
|
||||
*
|
||||
|
|
|
@ -76,7 +76,7 @@
|
|||
$ W( NMAX ), X( NMAX )
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAESM, CHKXER, DGELQ2, DGELQF, DGELQS, DORGL2,
|
||||
EXTERNAL ALAESM, CHKXER, DGELQ2, DGELQF, DORGL2,
|
||||
$ DORGLQ, DORML2, DORMLQ
|
||||
* ..
|
||||
* .. Scalars in Common ..
|
||||
|
@ -140,31 +140,6 @@
|
|||
CALL DGELQ2( 2, 1, A, 1, B, W, INFO )
|
||||
CALL CHKXER( 'DGELQ2', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* DGELQS
|
||||
*
|
||||
SRNAMT = 'DGELQS'
|
||||
INFOT = 1
|
||||
CALL DGELQS( -1, 0, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'DGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL DGELQS( 0, -1, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'DGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL DGELQS( 2, 1, 0, A, 2, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'DGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 3
|
||||
CALL DGELQS( 0, 0, -1, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'DGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 5
|
||||
CALL DGELQS( 2, 2, 0, A, 1, X, B, 2, W, 1, INFO )
|
||||
CALL CHKXER( 'DGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 8
|
||||
CALL DGELQS( 1, 2, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'DGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 10
|
||||
CALL DGELQS( 1, 1, 2, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'DGELQS', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* DORGLQ
|
||||
*
|
||||
SRNAMT = 'DORGLQ'
|
||||
|
|
|
@ -77,7 +77,7 @@
|
|||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAESM, CHKXER, DGEQR2, DGEQR2P, DGEQRF,
|
||||
$ DGEQRFP, DGEQRS, DORG2R, DORGQR, DORM2R,
|
||||
$ DGEQRFP, DORG2R, DORGQR, DORM2R,
|
||||
$ DORMQR
|
||||
* ..
|
||||
* .. Scalars in Common ..
|
||||
|
@ -170,31 +170,6 @@
|
|||
CALL DGEQR2P( 2, 1, A, 1, B, W, INFO )
|
||||
CALL CHKXER( 'DGEQR2P', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* DGEQRS
|
||||
*
|
||||
SRNAMT = 'DGEQRS'
|
||||
INFOT = 1
|
||||
CALL DGEQRS( -1, 0, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'DGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL DGEQRS( 0, -1, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'DGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL DGEQRS( 1, 2, 0, A, 2, X, B, 2, W, 1, INFO )
|
||||
CALL CHKXER( 'DGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 3
|
||||
CALL DGEQRS( 0, 0, -1, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'DGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 5
|
||||
CALL DGEQRS( 2, 1, 0, A, 1, X, B, 2, W, 1, INFO )
|
||||
CALL CHKXER( 'DGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 8
|
||||
CALL DGEQRS( 2, 1, 0, A, 2, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'DGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 10
|
||||
CALL DGEQRS( 1, 1, 2, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'DGEQRS', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* DORGQR
|
||||
*
|
||||
SRNAMT = 'DORGQR'
|
||||
|
|
|
@ -235,7 +235,7 @@
|
|||
REAL RESULT( NTESTS )
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, SERRLQ, SGELQS, SGET02,
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, SERRLQ, SGET02,
|
||||
$ SLACPY, SLARHS, SLATB4, SLATMS, SLQT01, SLQT02,
|
||||
$ SLQT03, XLAENV
|
||||
* ..
|
||||
|
@ -370,7 +370,7 @@
|
|||
$ WORK, LWORK, RWORK, RESULT( 3 ) )
|
||||
NT = NT + 4
|
||||
*
|
||||
* If M>=N and K=N, call SGELQS to solve a system
|
||||
* If M<=N and K=M, call SGELS to solve a system
|
||||
* with NRHS right hand sides and compute the
|
||||
* residual.
|
||||
*
|
||||
|
@ -387,14 +387,20 @@
|
|||
*
|
||||
CALL SLACPY( 'Full', M, NRHS, B, LDA, X,
|
||||
$ LDA )
|
||||
SRNAMT = 'SGELQS'
|
||||
CALL SGELQS( M, N, NRHS, AF, LDA, TAU, X,
|
||||
$ LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from SGELQS.
|
||||
* Reset AF to the original matrix. SGELS
|
||||
* factors the matrix before solving the system.
|
||||
*
|
||||
CALL SLACPY( 'Full', M, N, A, LDA, AF, LDA )
|
||||
*
|
||||
SRNAMT = 'SGELS'
|
||||
CALL SGELS( 'No transpose', M, N, NRHS, AF,
|
||||
$ LDA, X, LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from SGELS.
|
||||
*
|
||||
IF( INFO.NE.0 )
|
||||
$ CALL ALAERH( PATH, 'SGELQS', INFO, 0, ' ',
|
||||
$ CALL ALAERH( PATH, 'SGELS', INFO, 0, 'N',
|
||||
$ M, N, NRHS, -1, NB, IMAT,
|
||||
$ NFAIL, NERRS, NOUT )
|
||||
*
|
||||
|
|
|
@ -244,7 +244,7 @@
|
|||
EXTERNAL SGENND
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, SERRQR, SGEQRS, SGET02,
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, SERRQR, SGELS, SGET02,
|
||||
$ SLACPY, SLARHS, SLATB4, SLATMS, SQRT01,
|
||||
$ SQRT01P, SQRT02, SQRT03, XLAENV
|
||||
* ..
|
||||
|
@ -388,7 +388,7 @@
|
|||
$ WORK, LWORK, RWORK, RESULT( 3 ) )
|
||||
NT = NT + 4
|
||||
*
|
||||
* If M>=N and K=N, call SGEQRS to solve a system
|
||||
* If M>=N and K=N, call SGELS to solve a system
|
||||
* with NRHS right hand sides and compute the
|
||||
* residual.
|
||||
*
|
||||
|
@ -405,14 +405,20 @@
|
|||
*
|
||||
CALL SLACPY( 'Full', M, NRHS, B, LDA, X,
|
||||
$ LDA )
|
||||
SRNAMT = 'SGEQRS'
|
||||
CALL SGEQRS( M, N, NRHS, AF, LDA, TAU, X,
|
||||
$ LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from SGEQRS.
|
||||
* Reset AF to the original matrix. SGELS
|
||||
* factors the matrix before solving the system.
|
||||
*
|
||||
CALL SLACPY( 'Full', M, N, A, LDA, AF, LDA )
|
||||
*
|
||||
SRNAMT = 'SGELS'
|
||||
CALL SGELS( 'No transpose', M, N, NRHS, AF,
|
||||
$ LDA, X, LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from SGELS.
|
||||
*
|
||||
IF( INFO.NE.0 )
|
||||
$ CALL ALAERH( PATH, 'SGEQRS', INFO, 0, ' ',
|
||||
$ CALL ALAERH( PATH, 'SGELS', INFO, 0, 'N',
|
||||
$ M, N, NRHS, -1, NB, IMAT,
|
||||
$ NFAIL, NERRS, NOUT )
|
||||
*
|
||||
|
|
|
@ -76,7 +76,7 @@
|
|||
$ W( NMAX ), X( NMAX )
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAESM, CHKXER, SGELQ2, SGELQF, SGELQS, SORGL2,
|
||||
EXTERNAL ALAESM, CHKXER, SGELQ2, SGELQF, SORGL2,
|
||||
$ SORGLQ, SORML2, SORMLQ
|
||||
* ..
|
||||
* .. Scalars in Common ..
|
||||
|
@ -140,31 +140,6 @@
|
|||
CALL SGELQ2( 2, 1, A, 1, B, W, INFO )
|
||||
CALL CHKXER( 'SGELQ2', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* SGELQS
|
||||
*
|
||||
SRNAMT = 'SGELQS'
|
||||
INFOT = 1
|
||||
CALL SGELQS( -1, 0, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'SGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL SGELQS( 0, -1, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'SGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL SGELQS( 2, 1, 0, A, 2, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'SGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 3
|
||||
CALL SGELQS( 0, 0, -1, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'SGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 5
|
||||
CALL SGELQS( 2, 2, 0, A, 1, X, B, 2, W, 1, INFO )
|
||||
CALL CHKXER( 'SGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 8
|
||||
CALL SGELQS( 1, 2, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'SGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 10
|
||||
CALL SGELQS( 1, 1, 2, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'SGELQS', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* SORGLQ
|
||||
*
|
||||
SRNAMT = 'SORGLQ'
|
||||
|
|
|
@ -77,7 +77,7 @@
|
|||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAESM, CHKXER, SGEQR2, SGEQR2P, SGEQRF,
|
||||
$ SGEQRFP, SGEQRS, SORG2R, SORGQR, SORM2R,
|
||||
$ SGEQRFP, SORG2R, SORGQR, SORM2R,
|
||||
$ SORMQR
|
||||
* ..
|
||||
* .. Scalars in Common ..
|
||||
|
@ -170,31 +170,6 @@
|
|||
CALL SGEQR2P( 2, 1, A, 1, B, W, INFO )
|
||||
CALL CHKXER( 'SGEQR2P', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* SGEQRS
|
||||
*
|
||||
SRNAMT = 'SGEQRS'
|
||||
INFOT = 1
|
||||
CALL SGEQRS( -1, 0, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'SGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL SGEQRS( 0, -1, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'SGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL SGEQRS( 1, 2, 0, A, 2, X, B, 2, W, 1, INFO )
|
||||
CALL CHKXER( 'SGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 3
|
||||
CALL SGEQRS( 0, 0, -1, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'SGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 5
|
||||
CALL SGEQRS( 2, 1, 0, A, 1, X, B, 2, W, 1, INFO )
|
||||
CALL CHKXER( 'SGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 8
|
||||
CALL SGEQRS( 2, 1, 0, A, 2, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'SGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 10
|
||||
CALL SGEQRS( 1, 1, 2, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'SGEQRS', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* SORGQR
|
||||
*
|
||||
SRNAMT = 'SORGQR'
|
||||
|
|
|
@ -235,7 +235,7 @@
|
|||
DOUBLE PRECISION RESULT( NTESTS )
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRLQ, ZGELQS,
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRLQ, ZGELS,
|
||||
$ ZGET02, ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZLQT01,
|
||||
$ ZLQT02, ZLQT03
|
||||
* ..
|
||||
|
@ -370,7 +370,7 @@
|
|||
$ WORK, LWORK, RWORK, RESULT( 3 ) )
|
||||
NT = NT + 4
|
||||
*
|
||||
* If M>=N and K=N, call ZGELQS to solve a system
|
||||
* If M<=N and K=M, call ZGELS to solve a system
|
||||
* with NRHS right hand sides and compute the
|
||||
* residual.
|
||||
*
|
||||
|
@ -387,14 +387,20 @@
|
|||
*
|
||||
CALL ZLACPY( 'Full', M, NRHS, B, LDA, X,
|
||||
$ LDA )
|
||||
SRNAMT = 'ZGELQS'
|
||||
CALL ZGELQS( M, N, NRHS, AF, LDA, TAU, X,
|
||||
$ LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from ZGELQS.
|
||||
* Reset AF to the original matrix. ZGELS
|
||||
* factors the matrix before solving the system.
|
||||
*
|
||||
CALL ZLACPY( 'Full', M, N, A, LDA, AF, LDA )
|
||||
*
|
||||
SRNAMT = 'ZGELS'
|
||||
CALL ZGELS( 'No transpose', M, N, NRHS, AF,
|
||||
$ LDA, X, LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from ZGELS.
|
||||
*
|
||||
IF( INFO.NE.0 )
|
||||
$ CALL ALAERH( PATH, 'ZGELQS', INFO, 0, ' ',
|
||||
$ CALL ALAERH( PATH, 'ZGELS', INFO, 0, 'N',
|
||||
$ M, N, NRHS, -1, NB, IMAT,
|
||||
$ NFAIL, NERRS, NOUT )
|
||||
*
|
||||
|
|
|
@ -244,7 +244,7 @@
|
|||
EXTERNAL ZGENND
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRQR, ZGEQRS,
|
||||
EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRQR, ZGELS,
|
||||
$ ZGET02, ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZQRT01,
|
||||
$ ZQRT01P, ZQRT02, ZQRT03
|
||||
* ..
|
||||
|
@ -388,7 +388,7 @@
|
|||
$ WORK, LWORK, RWORK, RESULT( 3 ) )
|
||||
NT = NT + 4
|
||||
*
|
||||
* If M>=N and K=N, call ZGEQRS to solve a system
|
||||
* If M>=N and K=N, call ZGELS to solve a system
|
||||
* with NRHS right hand sides and compute the
|
||||
* residual.
|
||||
*
|
||||
|
@ -405,14 +405,20 @@
|
|||
*
|
||||
CALL ZLACPY( 'Full', M, NRHS, B, LDA, X,
|
||||
$ LDA )
|
||||
SRNAMT = 'ZGEQRS'
|
||||
CALL ZGEQRS( M, N, NRHS, AF, LDA, TAU, X,
|
||||
$ LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from ZGEQRS.
|
||||
* Reset AF to the original matrix. ZGELS
|
||||
* factors the matrix before solving the system.
|
||||
*
|
||||
CALL ZLACPY( 'Full', M, N, A, LDA, AF, LDA )
|
||||
*
|
||||
SRNAMT = 'ZGELS'
|
||||
CALL ZGELS( 'No transpose', M, N, NRHS, AF,
|
||||
$ LDA, X, LDA, WORK, LWORK, INFO )
|
||||
*
|
||||
* Check error code from ZGELS.
|
||||
*
|
||||
IF( INFO.NE.0 )
|
||||
$ CALL ALAERH( PATH, 'ZGEQRS', INFO, 0, ' ',
|
||||
$ CALL ALAERH( PATH, 'ZGELS', INFO, 0, 'N',
|
||||
$ M, N, NRHS, -1, NB, IMAT,
|
||||
$ NFAIL, NERRS, NOUT )
|
||||
*
|
||||
|
|
|
@ -76,7 +76,7 @@
|
|||
$ W( NMAX ), X( NMAX )
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAESM, CHKXER, ZGELQ2, ZGELQF, ZGELQS, ZUNGL2,
|
||||
EXTERNAL ALAESM, CHKXER, ZGELQ2, ZGELQF, ZUNGL2,
|
||||
$ ZUNGLQ, ZUNML2, ZUNMLQ
|
||||
* ..
|
||||
* .. Scalars in Common ..
|
||||
|
@ -142,31 +142,6 @@
|
|||
CALL ZGELQ2( 2, 1, A, 1, B, W, INFO )
|
||||
CALL CHKXER( 'ZGELQ2', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* ZGELQS
|
||||
*
|
||||
SRNAMT = 'ZGELQS'
|
||||
INFOT = 1
|
||||
CALL ZGELQS( -1, 0, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL ZGELQS( 0, -1, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL ZGELQS( 2, 1, 0, A, 2, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 3
|
||||
CALL ZGELQS( 0, 0, -1, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 5
|
||||
CALL ZGELQS( 2, 2, 0, A, 1, X, B, 2, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 8
|
||||
CALL ZGELQS( 1, 2, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGELQS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 10
|
||||
CALL ZGELQS( 1, 1, 2, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGELQS', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* ZUNGLQ
|
||||
*
|
||||
SRNAMT = 'ZUNGLQ'
|
||||
|
|
|
@ -77,7 +77,7 @@
|
|||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL ALAESM, CHKXER, ZGEQR2, ZGEQR2P, ZGEQRF,
|
||||
$ ZGEQRFP, ZGEQRS, ZUNG2R, ZUNGQR, ZUNM2R,
|
||||
$ ZGEQRFP, ZUNG2R, ZUNGQR, ZUNM2R,
|
||||
$ ZUNMQR
|
||||
* ..
|
||||
* .. Scalars in Common ..
|
||||
|
@ -172,31 +172,6 @@
|
|||
CALL ZGEQR2P( 2, 1, A, 1, B, W, INFO )
|
||||
CALL CHKXER( 'ZGEQR2P', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* ZGEQRS
|
||||
*
|
||||
SRNAMT = 'ZGEQRS'
|
||||
INFOT = 1
|
||||
CALL ZGEQRS( -1, 0, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL ZGEQRS( 0, -1, 0, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 2
|
||||
CALL ZGEQRS( 1, 2, 0, A, 2, X, B, 2, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 3
|
||||
CALL ZGEQRS( 0, 0, -1, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 5
|
||||
CALL ZGEQRS( 2, 1, 0, A, 1, X, B, 2, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 8
|
||||
CALL ZGEQRS( 2, 1, 0, A, 2, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGEQRS', INFOT, NOUT, LERR, OK )
|
||||
INFOT = 10
|
||||
CALL ZGEQRS( 1, 1, 2, A, 1, X, B, 1, W, 1, INFO )
|
||||
CALL CHKXER( 'ZGEQRS', INFOT, NOUT, LERR, OK )
|
||||
*
|
||||
* ZUNGQR
|
||||
*
|
||||
SRNAMT = 'ZUNGQR'
|
||||
|
|
Loading…
Reference in New Issue