floraachy
/
OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #2797 from martin-frbg/relafixes1 

ReLAPACK fixes
This commit is contained in:
Martin Kroeker 2020-09-01 16:04:03 +02:00 committed by GitHub
parent f42e84d46c de63675717
commit 4074770d00
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
33 changed files with 163 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
relapack/src/cgbtrf.c
View File

@ -36,6 +36,7 @@ void RELAPACK_cgbtrf(
return;
}
if (*m == 0 || *n == 0) return;
// Constant
const float ZERO[] = { 0., 0. };
@ -56,10 +57,10 @@ void RELAPACK_cgbtrf(
// Allocate work space
const blasint n1 = CREC_SPLIT(*n);
const blasint mWorkl = (kv > n1) ? MAX(1, *m - *kl) : kv;
const blasint nWorkl = (kv > n1) ? n1 : kv;
const blasint mWorku = (*kl > n1) ? n1 : *kl;
const blasint nWorku = (*kl > n1) ? MAX(0, *n - *kl) : *kl;
const blasint mWorkl = abs ( (kv > n1) ? MAX(1, *m - *kl) : kv);
const blasint nWorkl = abs ( (kv > n1) ? n1 : kv);
const blasint mWorku = abs ((*kl > n1) ? n1 : *kl);
const blasint nWorku = abs ((*kl > n1) ? MAX(0, *n - *kl) : *kl);
float *Workl = malloc(mWorkl * nWorkl * 2 * sizeof(float));
float *Worku = malloc(mWorku * nWorku * 2 * sizeof(float));
LAPACK(claset)("L", &mWorkl, &nWorkl, ZERO, ZERO, Workl, &mWorkl);
@ -82,7 +83,7 @@ static void RELAPACK_cgbtrf_rec(
blasint *info
) {
if (*n <= MAX(CROSSOVER_CGBTRF, 1)) {
if (*n <= MAX(CROSSOVER_CGBTRF, 1)|| *n > *kl || *ldAb == 1) {
// Unblocked
LAPACK(cgbtf2)(m, n, kl, ku, Ab, ldAb, ipiv, info);
return;

9
relapack/src/cgetrf.c
View File

@ -30,6 +30,8 @@ void RELAPACK_cgetrf(
return;
}
if (*m == 0 || *n == 0) return;
const blasint sn = MIN(*m, *n);
RELAPACK_cgetrf_rec(m, &sn, A, ldA, ipiv, info);
@ -62,9 +64,11 @@ static void RELAPACK_cgetrf_rec(
blasint *info
) {
if (*n <= MAX(CROSSOVER_CGETRF, 1)) {
if (*m == 0 || *n == 0) return;
if ( *n <= MAX(CROSSOVER_CGETRF, 1)) {
// Unblocked
LAPACK(cgetf2)(m, n, A, ldA, ipiv, info);
LAPACK(cgetrf2)(m, n, A, ldA, ipiv, info);
return;
}
@ -96,6 +100,7 @@ static void RELAPACK_cgetrf_rec(
// recursion(A_L, ipiv_T)
RELAPACK_cgetrf_rec(m, &n1, A_L, ldA, ipiv_T, info);
if (*info) return;
// apply pivots to A_R
LAPACK(claswp)(&n2, A_R, ldA, iONE, &n1, ipiv_T, iONE);

2
relapack/src/chegst.c
View File

@ -40,6 +40,8 @@ void RELAPACK_chegst(
return;
}
if (*n == 0) return;
// Clean char * arguments
const char cleanuplo = lower ? 'L' : 'U';

4
relapack/src/chetrf_rook.c
View File

@ -36,7 +36,7 @@ void RELAPACK_chetrf_rook(
*info = -2;
else if (*ldA < MAX(1, *n))
*info = -4;
else if (*lWork < minlWork && *lWork != -1)
else if ((*lWork < 1 || *lWork < minlWork) && *lWork != -1)
*info = -7;
else if (*lWork == -1) {
// Work size query
@ -56,7 +56,7 @@ void RELAPACK_chetrf_rook(
if (*info) {
const blasint minfo = -*info;
LAPACK(xerbla)("CHETRF", &minfo, strlen("CHETRF"));
LAPACK(xerbla)("CHETRF_ROOK", &minfo, strlen("CHETRF_ROOK"));
return;
}

2
relapack/src/clauum.c
View File

@ -32,6 +32,8 @@ void RELAPACK_clauum(
return;
}
if (*n == 0) return;
// Clean char * arguments
const char cleanuplo = lower ? 'L' : 'U';

10
relapack/src/cpbtrf.c
View File

@ -35,6 +35,8 @@ void RELAPACK_cpbtrf(
return;
}
if (*n == 0) return;
// Clean char * arguments
const char cleanuplo = lower ? 'L' : 'U';
@ -43,8 +45,8 @@ void RELAPACK_cpbtrf(
// Allocate work space
const blasint n1 = CREC_SPLIT(*n);
const blasint mWork = (*kd > n1) ? (lower ? *n - *kd : n1) : *kd;
const blasint nWork = (*kd > n1) ? (lower ? n1 : *n - *kd) : *kd;
const blasint mWork = abs((*kd > n1) ? (lower ? *n - *kd : n1) : *kd);
const blasint nWork = abs((*kd > n1) ? (lower ? n1 : *n - *kd) : *kd);
float *Work = malloc(mWork * nWork * 2 * sizeof(float));
LAPACK(claset)(uplo, &mWork, &nWork, ZERO, ZERO, Work, &mWork);
@ -64,7 +66,7 @@ static void RELAPACK_cpbtrf_rec(
blasint *info
){
if (*n <= MAX(CROSSOVER_CPBTRF, 1)) {
if (*n <= MAX(CROSSOVER_CPBTRF, 1) || *ldAb==1) {
// Unblocked
LAPACK(cpbtf2)(uplo, n, kd, Ab, ldAb, info);
return;
@ -148,7 +150,7 @@ static void RELAPACK_cpbtrf_rec(
}
// recursion(A_BR)
if (*kd > n1)
if (*kd > n1 && ldA != 0)
RELAPACK_cpotrf(uplo, &n2, A_BR, ldA, info);
else
RELAPACK_cpbtrf_rec(uplo, &n2, kd, Ab_BR, ldAb, Work, ldWork, info);

3
relapack/src/cpotrf.c
View File

@ -32,6 +32,8 @@ void RELAPACK_cpotrf(
return;
}
if (*n == 0) return;
// Clean char * arguments
const char cleanuplo = lower ? 'L' : 'U';
@ -46,6 +48,7 @@ static void RELAPACK_cpotrf_rec(
float *A, const blasint *ldA,
blasint *info
){
if (*n == 0) return;
if (*n <= MAX(CROSSOVER_CPOTRF, 1)) {
// Unblocked

3
relapack/src/csytrf.c
View File

@ -36,7 +36,7 @@ void RELAPACK_csytrf(
*info = -2;
else if (*ldA < MAX(1, *n))
*info = -4;
else if (*lWork < minlWork && *lWork != -1)
else if ((*lWork < 1 || *lWork < minlWork) && *lWork != -1)
*info = -7;
else if (*lWork == -1) {
// Work size query
@ -67,6 +67,7 @@ void RELAPACK_csytrf(
blasint nout;
// Recursive kernel
if (*n != 0)
RELAPACK_csytrf_rec(&cleanuplo, n, n, &nout, A, ldA, ipiv, cleanWork, n, info);
#if XSYTRF_ALLOW_MALLOC

4
relapack/src/csytrf_rook.c
View File

@ -36,7 +36,7 @@ void RELAPACK_csytrf_rook(
*info = -2;
else if (*ldA < MAX(1, *n))
*info = -4;
else if (*lWork < minlWork && *lWork != -1)
else if ((*lWork < 1 || *lWork < minlWork) && *lWork != -1)
*info = -7;
else if (*lWork == -1) {
// Work size query
@ -56,7 +56,7 @@ void RELAPACK_csytrf_rook(
if (*info) {
const blasint minfo = -*info;
LAPACK(xerbla)("CSYTRF", &minfo, strlen("CSYTRF"));
LAPACK(xerbla)("CSYTRF_ROOK", &minfo, strlen("CSYTRF_ROOK"));
return;
}

7
relapack/src/ctgsyl.c
View File

@ -68,6 +68,13 @@ void RELAPACK_ctgsyl(
return;
}
if ( *m == 0 || *n == 0) {
*scale = 1.;
if (notran && (*ijob != 0))
*dif = 0.;
return;
}
// Clean char * arguments
const char cleantrans = notran ? 'N' : 'C';

5
relapack/src/ctrsyl.c
View File

@ -47,6 +47,11 @@ void RELAPACK_ctrsyl(
return;
}
if (*m == 0 || *n == 0) {
*scale = 1.;
return;
}
// Clean char * arguments
const char cleantranA = notransA ? 'N' : 'C';
const char cleantranB = notransB ? 'N' : 'C';

2
relapack/src/ctrtri.c
View File

@ -36,6 +36,8 @@ void RELAPACK_ctrtri(
return;
}
if (*n == 0) return;
// Clean char * arguments
const char cleanuplo = lower ? 'L' : 'U';
const char cleandiag = nounit ? 'N' : 'U';

5
relapack/src/dgbtrf.c
View File

@ -36,6 +36,8 @@ void RELAPACK_dgbtrf(
return;
}
if (*m == 0 || *n == 0) return;
// Constant
const double ZERO[] = { 0. };
@ -83,7 +85,7 @@ static void RELAPACK_dgbtrf_rec(
blasint *info
) {
if (*n <= MAX(CROSSOVER_DGBTRF, 1)) {
if (*n <= MAX(CROSSOVER_DGBTRF, 1) || *n > *kl || *ldAb == 1) {
// Unblocked
LAPACK(dgbtf2)(m, n, kl, ku, Ab, ldAb, ipiv, info);
return;
@ -195,6 +197,7 @@ static void RELAPACK_dgbtrf_rec(
// Worku = A_TRr
LAPACK(dlacpy)("L", &m1, &n22, A_TRr, ldA, Worku, ldWorku);
// Worku = A_TL \ Worku
if (ldWorku <= 0) return;
BLAS(dtrsm)("L", "L", "N", "U", &m1, &n22, ONE, A_TL, ldA, Worku, ldWorku);
// A_TRr = Worku
LAPACK(dlacpy)("L", &m1, &n22, Worku, ldWorku, A_TRr, ldA);

12
relapack/src/dgetrf.c
View File

@ -29,15 +29,16 @@ void RELAPACK_dgetrf(
return;
}
const blasint sn = MIN(*m, *n);
if (*m == 0 || *n == 0) return;
const blasint sn = MIN(*m, *n);
RELAPACK_dgetrf_rec(m, &sn, A, ldA, ipiv, info);
// Right remainder
if (*m < *n) {
// Constants
const double ONE[] = { 1. };
const blasint iONE[] = { 1. };
const blasint iONE[] = { 1 };
// Splitting
const blasint rn = *n - *m;
@ -60,13 +61,11 @@ static void RELAPACK_dgetrf_rec(
double *A, const blasint *ldA, blasint *ipiv,
blasint *info
) {
if (*n <= MAX(CROSSOVER_DGETRF, 1)) {
if ( *n <= MAX(CROSSOVER_DGETRF, 1)) {
// Unblocked
LAPACK(dgetf2)(m, n, A, ldA, ipiv, info);
LAPACK(dgetrf2)(m, n, A, ldA, ipiv, info);
return;
}
// Constants
const double ONE[] = { 1. };
const double MONE[] = { -1. };
@ -95,6 +94,7 @@ static void RELAPACK_dgetrf_rec(
// recursion(A_L, ipiv_T)
RELAPACK_dgetrf_rec(m, &n1, A_L, ldA, ipiv_T, info);
if (*info) return;
// apply pivots to A_R
LAPACK(dlaswp)(&n2, A_R, ldA, iONE, &n1, ipiv_T, iONE);

10
relapack/src/dpbtrf.c
View File

@ -35,6 +35,8 @@ void RELAPACK_dpbtrf(
return;
}
if (*n == 0) return;
// Clean char * arguments
const char cleanuplo = lower ? 'L' : 'U';
@ -43,8 +45,8 @@ void RELAPACK_dpbtrf(
// Allocate work space
const blasint n1 = DREC_SPLIT(*n);
const blasint mWork = (*kd > n1) ? (lower ? *n - *kd : n1) : *kd;
const blasint nWork = (*kd > n1) ? (lower ? n1 : *n - *kd) : *kd;
const blasint mWork = abs((*kd > n1) ? (lower ? *n - *kd : n1) : *kd);
const blasint nWork = abs((*kd > n1) ? (lower ? n1 : *n - *kd) : *kd);
double *Work = malloc(mWork * nWork * sizeof(double));
LAPACK(dlaset)(uplo, &mWork, &nWork, ZERO, ZERO, Work, &mWork);
@ -64,7 +66,7 @@ static void RELAPACK_dpbtrf_rec(
blasint *info
){
if (*n <= MAX(CROSSOVER_DPBTRF, 1)) {
if (*n <= MAX(CROSSOVER_DPBTRF, 1) || *ldAb == 1) {
// Unblocked
LAPACK(dpbtf2)(uplo, n, kd, Ab, ldAb, info);
return;
@ -148,7 +150,7 @@ static void RELAPACK_dpbtrf_rec(
}
// recursion(A_BR)
if (*kd > n1)
if (*kd > n1 && ldA != 0)
RELAPACK_dpotrf(uplo, &n2, A_BR, ldA, info);
else
RELAPACK_dpbtrf_rec(uplo, &n2, kd, Ab_BR, ldAb, Work, ldWork, info);

3
relapack/src/dsytrf.c
View File

@ -36,7 +36,7 @@ void RELAPACK_dsytrf(
*info = -2;
else if (*ldA < MAX(1, *n))
*info = -4;
else if (*lWork < minlWork && *lWork != -1)
else if ((*lWork < 1 || *lWork < minlWork) && *lWork != -1)
*info = -7;
else if (*lWork == -1) {
// Work size query
@ -67,6 +67,7 @@ void RELAPACK_dsytrf(
blasint nout;
// Recursive kernel
if (*n != 0)
RELAPACK_dsytrf_rec(&cleanuplo, n, n, &nout, A, ldA, ipiv, cleanWork, n, info);
#if XSYTRF_ALLOW_MALLOC

4
relapack/src/dsytrf_rook.c
View File

@ -36,7 +36,7 @@ void RELAPACK_dsytrf_rook(
*info = -2;
else if (*ldA < MAX(1, *n))
*info = -4;
else if (*lWork < minlWork && *lWork != -1)
else if ((*lWork <1 || *lWork < minlWork) && *lWork != -1)
*info = -7;
else if (*lWork == -1) {
// Work size query
@ -56,7 +56,7 @@ void RELAPACK_dsytrf_rook(
if (*info) {
const blasint minfo = -*info;
LAPACK(xerbla)("DSYTRF", &minfo, strlen("DSYTRF"));
LAPACK(xerbla)("DSYTRF_ROOK", &minfo, strlen("DSYTRF_ROOK"));
return;
}

5
relapack/src/dtrsyl.c
View File

@ -49,6 +49,11 @@ void RELAPACK_dtrsyl(
return;
}
if (*m == 0 || *n == 0) {
*scale = 1.;
return;
}
// Clean char * arguments
const char cleantranA = notransA ? 'N' : (transA ? 'T' : 'C');
const char cleantranB = notransB ? 'N' : (transB ? 'T' : 'C');

7
relapack/src/lapack.h
View File

@ -4,6 +4,13 @@
extern blasint LAPACK(lsame)(const char *, const char *);
extern blasint LAPACK(xerbla)(const char *, const blasint *, int);
extern const blasint LAPACK(ilaenv)(const blasint *, const char*, const char*, const blasint* , int , int, int );
extern void LAPACK(sgetrf2)(const blasint *, const blasint *, float *, const blasint *, blasint *, blasint *);
extern void LAPACK(dgetrf2)(const blasint *, const blasint *, double *, const blasint *, blasint *, blasint *);
extern void LAPACK(cgetrf2)(const blasint *, const blasint *, float *, const blasint *, blasint *, blasint *);
extern void LAPACK(zgetrf2)(const blasint *, const blasint *, double *, const blasint *, blasint *, blasint *);
extern void LAPACK(slaswp)(const blasint *, float *, const blasint *, const blasint *, const blasint *, const blasint *, const blasint *);
extern void LAPACK(dlaswp)(const blasint *, double *, const blasint *, const blasint *, const blasint *, const blasint *, const blasint *);
extern void LAPACK(claswp)(const blasint *, float *, const blasint *, const blasint *, const blasint *, const blasint *, const blasint *);

16
relapack/src/sgbtrf.c
View File

@ -35,6 +35,13 @@ void RELAPACK_sgbtrf(
return;
}
if (*m == 0 || *n == 0) return;
if (*ldAb == 1) {
LAPACK(sgbtf2)(m, n, kl, ku, Ab, ldAb, ipiv, info);
return;
}
// Constant
const float ZERO[] = { 0. };
@ -82,8 +89,9 @@ static void RELAPACK_sgbtrf_rec(
blasint *info
) {
if (*m == 0 || *n == 0) return;
if (*n <= MAX(CROSSOVER_SGBTRF, 1)) {
if ( *n <= MAX(CROSSOVER_SGBTRF, 1) || *n > *kl || *ldAb == 1) {
// Unblocked
LAPACK(sgbtf2)(m, n, kl, ku, Ab, ldAb, ipiv, info);
return;
@ -160,7 +168,7 @@ static void RELAPACK_sgbtrf_rec(
// recursion(Ab_L, ipiv_T)
RELAPACK_sgbtrf_rec(m, &n1, kl, ku, Ab_L, ldAb, ipiv_T, Workl, ldWorkl, Worku, ldWorku, info);
if (*info) return;
// Workl = A_BLb
LAPACK(slacpy)("U", &m22, &n1, A_BLb, ldA, Workl, ldWorkl);
@ -222,8 +230,8 @@ static void RELAPACK_sgbtrf_rec(
// recursion(Ab_BR, ipiv_B)
//cause of infinite recursion here ?
// RELAPACK_sgbtrf_rec(&m2, &n2, kl, ku, Ab_BR, ldAb, ipiv_B, Workl, ldWorkl, Worku, ldWorku, info);
LAPACK(sgbtf2)(&m2, &n2, kl, ku, Ab_BR, ldAb, ipiv_B, info);
RELAPACK_sgbtrf_rec(&m2, &n2, kl, ku, Ab_BR, ldAb, ipiv_B, Workl, ldWorkl, Worku, ldWorku, info);
// LAPACK(sgbtf2)(&m2, &n2, kl, ku, Ab_BR, ldAb, ipiv_B, info);
if (*info)
*info += n1;
// shift pivots

15
relapack/src/sgetrf.c
View File

@ -14,7 +14,6 @@ void RELAPACK_sgetrf(
float *A, const blasint *ldA, blasint *ipiv,
blasint *info
) {
// Check arguments
*info = 0;
if (*m < 0)
@ -28,6 +27,9 @@ void RELAPACK_sgetrf(
LAPACK(xerbla)("SGETRF", &minfo, strlen("SGETRF"));
return;
}
if (*m == 0 || *n == 0) return;
const blasint sn = MIN(*m, *n);
RELAPACK_sgetrf_rec(m, &sn, A, ldA, ipiv, info);
@ -35,7 +37,7 @@ void RELAPACK_sgetrf(
if (*m < *n) {
// Constants
const float ONE[] = { 1. };
const blasint iONE[] = { 1. };
const blasint iONE[] = { 1 };
// Splitting
const blasint rn = *n - *m;
@ -58,9 +60,12 @@ static void RELAPACK_sgetrf_rec(
float *A, const blasint *ldA, blasint *ipiv,
blasint *info
) {
if (*n <= MAX(CROSSOVER_SGETRF, 1)) {
if (*m == 0 || *n == 0) return;
if ( *n <= MAX(CROSSOVER_SGETRF, 1)) {
// Unblocked
LAPACK(sgetf2)(m, n, A, ldA, ipiv, info);
LAPACK(sgetrf2)(m, n, A, ldA, ipiv, info);
return;
}
@ -91,6 +96,8 @@ static void RELAPACK_sgetrf_rec(
// recursion(A_L, ipiv_T)
RELAPACK_sgetrf_rec(m, &n1, A_L, ldA, ipiv_T, info);
if (*info)
return;
// apply pivots to A_R
LAPACK(slaswp)(&n2, A_R, ldA, iONE, &n1, ipiv_T, iONE);

13
relapack/src/spbtrf.c
View File

@ -35,6 +35,9 @@ void RELAPACK_spbtrf(
return;
}
if (*n == 0) return;
// Clean char * arguments
const char cleanuplo = lower ? 'L' : 'U';
@ -43,8 +46,8 @@ void RELAPACK_spbtrf(
// Allocate work space
const blasint n1 = SREC_SPLIT(*n);
const blasint mWork = (*kd > n1) ? (lower ? *n - *kd : n1) : *kd;
const blasint nWork = (*kd > n1) ? (lower ? n1 : *n - *kd) : *kd;
const blasint mWork = abs( (*kd > n1) ? (lower ? *n - *kd : n1) : *kd);
const blasint nWork = abs((*kd > n1) ? (lower ? n1 : *n - *kd) : *kd);
float *Work = malloc(mWork * nWork * sizeof(float));
LAPACK(slaset)(uplo, &mWork, &nWork, ZERO, ZERO, Work, &mWork);
@ -64,7 +67,9 @@ static void RELAPACK_spbtrf_rec(
blasint *info
){
if (*n <= MAX(CROSSOVER_SPBTRF, 1)) {
if (*n == 0 ) return;
if ( *n <= MAX(CROSSOVER_SPBTRF, 1) || *ldAb == 1) {
// Unblocked
LAPACK(spbtf2)(uplo, n, kd, Ab, ldAb, info);
return;
@ -148,7 +153,7 @@ static void RELAPACK_spbtrf_rec(
}
// recursion(A_BR)
if (*kd > n1)
if (*kd > n1 && ldA != 0)
RELAPACK_spotrf(uplo, &n2, A_BR, ldA, info);
else
RELAPACK_spbtrf_rec(uplo, &n2, kd, Ab_BR, ldAb, Work, ldWork, info);

3
relapack/src/ssytrf.c
View File

@ -35,7 +35,7 @@ void RELAPACK_ssytrf(
*info = -2;
else if (*ldA < MAX(1, *n))
*info = -4;
else if (*lWork < minlWork && *lWork != -1)
else if ((*lWork <1 || *lWork < minlWork) && *lWork != -1)
*info = -7;
else if (*lWork == -1) {
// Work size query
@ -66,6 +66,7 @@ void RELAPACK_ssytrf(
blasint nout;
// Recursive kernel
if (*n != 0)
RELAPACK_ssytrf_rec(&cleanuplo, n, n, &nout, A, ldA, ipiv, cleanWork, n, info);
#if XSYTRF_ALLOW_MALLOC

5
relapack/src/ssytrf_rook.c
View File

@ -36,7 +36,7 @@ void RELAPACK_ssytrf_rook(
*info = -2;
else if (*ldA < MAX(1, *n))
*info = -4;
else if (*lWork < minlWork && *lWork != -1)
else if ((*lWork < 1 ||*lWork < minlWork) && *lWork != -1)
*info = -7;
else if (*lWork == -1) {
// Work size query
@ -56,7 +56,7 @@ void RELAPACK_ssytrf_rook(
if (*info) {
const blasint minfo = -*info;
LAPACK(xerbla)("SSYTRF", &minfo, strlen("SSYTRF"));
LAPACK(xerbla)("SSYTRF_ROOK", &minfo, strlen("SSYTRF_ROOK"));
return;
}
@ -67,6 +67,7 @@ void RELAPACK_ssytrf_rook(
blasint nout;
// Recursive kernel
if (*n != 0)
RELAPACK_ssytrf_rook_rec(&cleanuplo, n, n, &nout, A, ldA, ipiv, cleanWork, n, info);
#if XSYTRF_ALLOW_MALLOC

5
relapack/src/strsyl.c
View File

@ -49,6 +49,11 @@ void RELAPACK_strsyl(
return;
}
if (*m == 0 || *n == 0) {
*scale = 1.;
return;
}
// Clean char * arguments
const char cleantranA = notransA ? 'N' : (transA ? 'T' : 'C');
const char cleantranB = notransB ? 'N' : (transB ? 'T' : 'C');

16
relapack/src/zgbtrf.c
View File

@ -36,6 +36,8 @@ void RELAPACK_zgbtrf(
return;
}
if (*m == 0 || *n == 0) return;
// Constant
const double ZERO[] = { 0., 0. };
@ -82,7 +84,7 @@ static void RELAPACK_zgbtrf_rec(
blasint *info
) {
if (*n <= MAX(CROSSOVER_ZGBTRF, 1)) {
if (*n <= MAX(CROSSOVER_ZGBTRF, 1) || *n > *kl || *ldAb == 1) {
// Unblocked
LAPACK(zgbtf2)(m, n, kl, ku, Ab, ldAb, ipiv, info);
return;
@ -92,6 +94,7 @@ static void RELAPACK_zgbtrf_rec(
const double ONE[] = { 1., 0. };
const double MONE[] = { -1., 0. };
const blasint iONE[] = { 1 };
const blasint min11 = -11;
// Loop iterators
blasint i, j;
@ -158,6 +161,7 @@ static void RELAPACK_zgbtrf_rec(
// recursion(Ab_L, ipiv_T)
RELAPACK_zgbtrf_rec(m, &n1, kl, ku, Ab_L, ldAb, ipiv_T, Workl, ldWorkl, Worku, ldWorku, info);
if (*info) return;
// Workl = A_BLb
LAPACK(zlacpy)("U", &m22, &n1, A_BLb, ldA, Workl, ldWorkl);
@ -193,11 +197,21 @@ static void RELAPACK_zgbtrf_rec(
}
// A_TRl = A_TL \ A_TRl
if (*ldA < MAX(1,m1)) {
LAPACK(xerbla)("ZGBTRF", &min11, strlen("ZGBTRF"));
return;
} else {
BLAS(ztrsm)("L", "L", "N", "U", &m1, &n21, ONE, A_TL, ldA, A_TRl, ldA);
}
// Worku = A_TRr
LAPACK(zlacpy)("L", &m1, &n22, A_TRr, ldA, Worku, ldWorku);
// Worku = A_TL \ Worku
if (*ldWorku < MAX(1,m1)) {
LAPACK(xerbla)("ZGBTRF", &min11, strlen("ZGBTRF"));
return;
} else {
BLAS(ztrsm)("L", "L", "N", "U", &m1, &n22, ONE, A_TL, ldA, Worku, ldWorku);
}
// A_TRr = Worku
LAPACK(zlacpy)("L", &m1, &n22, Worku, ldWorku, A_TRr, ldA);
// A_BRtl = A_BRtl - A_BLt * A_TRl

9
relapack/src/zgetrf.c
View File

@ -30,6 +30,7 @@ void RELAPACK_zgetrf(
return;
}
if (*m == 0 || *n == 0) return;
const blasint sn = MIN(*m, *n);
RELAPACK_zgetrf_rec(m, &sn, A, ldA, ipiv, info);
@ -62,9 +63,11 @@ static void RELAPACK_zgetrf_rec(
blasint *info
) {
if (*n <= MAX(CROSSOVER_ZGETRF, 1)) {
if (*m == 0 || *n == 0) return;
if ( *n <= MAX(CROSSOVER_ZGETRF, 1)) {
// Unblocked
LAPACK(zgetf2)(m, n, A, ldA, ipiv, info);
LAPACK(zgetrf2)(m, n, A, ldA, ipiv, info);
return;
}
@ -96,6 +99,8 @@ static void RELAPACK_zgetrf_rec(
// recursion(A_L, ipiv_T)
RELAPACK_zgetrf_rec(m, &n1, A_L, ldA, ipiv_T, info);
if (*info) return;
// apply pivots to A_R
LAPACK(zlaswp)(&n2, A_R, ldA, iONE, &n1, ipiv_T, iONE);

4
relapack/src/zhetrf_rook.c
View File

@ -36,7 +36,7 @@ void RELAPACK_zhetrf_rook(
*info = -2;
else if (*ldA < MAX(1, *n))
*info = -4;
else if (*lWork < minlWork && *lWork != -1)
else if ((*lWork < 1 || *lWork < minlWork) && *lWork != -1)
*info = -7;
else if (*lWork == -1) {
// Work size query
@ -56,7 +56,7 @@ void RELAPACK_zhetrf_rook(
if (*info) {
const blasint minfo = -*info;
LAPACK(xerbla)("ZHETRF", &minfo, strlen("ZHETRF"));
LAPACK(xerbla)("ZHETRF_ROOK", &minfo, strlen("ZHETRF_ROOK"));
return;
}

11
relapack/src/zpbtrf.c
View File

@ -35,6 +35,8 @@ void RELAPACK_zpbtrf(
return;
}
if (*n == 0) return;
// Clean char * arguments
const char cleanuplo = lower ? 'L' : 'U';
@ -43,9 +45,10 @@ void RELAPACK_zpbtrf(
// Allocate work space
const blasint n1 = ZREC_SPLIT(*n);
const blasint mWork = (*kd > n1) ? (lower ? *n - *kd : n1) : *kd;
const blasint nWork = (*kd > n1) ? (lower ? n1 : *n - *kd) : *kd;
const blasint mWork = abs((*kd > n1) ? (lower ? *n - *kd : n1) : *kd);
const blasint nWork = abs((*kd > n1) ? (lower ? n1 : *n - *kd) : *kd);
double *Work = malloc(mWork * nWork * 2 * sizeof(double));
LAPACK(zlaset)(uplo, &mWork, &nWork, ZERO, ZERO, Work, &mWork);
// Recursive kernel
@ -64,7 +67,7 @@ static void RELAPACK_zpbtrf_rec(
blasint *info
){
if (*n <= MAX(CROSSOVER_ZPBTRF, 1)) {
if (*n <= MAX(CROSSOVER_ZPBTRF, 1) || *ldAb == 1) {
// Unblocked
LAPACK(zpbtf2)(uplo, n, kd, Ab, ldAb, info);
return;
@ -148,7 +151,7 @@ static void RELAPACK_zpbtrf_rec(
}
// recursion(A_BR)
if (*kd > n1)
if (*kd > n1 && ldA != 0)
RELAPACK_zpotrf(uplo, &n2, A_BR, ldA, info);
else
RELAPACK_zpbtrf_rec(uplo, &n2, kd, Ab_BR, ldAb, Work, ldWork, info);

3
relapack/src/zsytrf.c
View File

@ -36,7 +36,7 @@ void RELAPACK_zsytrf(
*info = -2;
else if (*ldA < MAX(1, *n))
*info = -4;
else if (*lWork < minlWork && *lWork != -1)
else if ((*lWork < 1 || *lWork < minlWork) && *lWork != -1)
*info = -7;
else if (*lWork == -1) {
// Work size query
@ -67,6 +67,7 @@ void RELAPACK_zsytrf(
blasint nout;
// Recursive kernel
if (*n != 0)
RELAPACK_zsytrf_rec(&cleanuplo, n, n, &nout, A, ldA, ipiv, cleanWork, n, info);
#if XSYTRF_ALLOW_MALLOC

5
relapack/src/zsytrf_rook.c
View File

@ -36,7 +36,7 @@ void RELAPACK_zsytrf_rook(
*info = -2;
else if (*ldA < MAX(1, *n))
*info = -4;
else if (*lWork < minlWork && *lWork != -1)
else if ((*lWork < 1 || *lWork < minlWork) && *lWork != -1)
*info = -7;
else if (*lWork == -1) {
// Work size query
@ -56,7 +56,7 @@ void RELAPACK_zsytrf_rook(
if (*info) {
const blasint minfo = -*info;
LAPACK(xerbla)("ZSYTRF", &minfo, strlen("ZSYTRF"));
LAPACK(xerbla)("ZSYTRF_ROOK", &minfo, strlen("ZSYTRF_ROOK"));
return;
}
@ -67,6 +67,7 @@ void RELAPACK_zsytrf_rook(
blasint nout;
// Recursive kernel
if (*n != 0)
RELAPACK_zsytrf_rook_rec(&cleanuplo, n, n, &nout, A, ldA, ipiv, cleanWork, n, info);
#if XSYTRF_ALLOW_MALLOC

5
relapack/src/ztrsyl.c
View File

@ -47,6 +47,11 @@ void RELAPACK_ztrsyl(
return;
}
if (*m == 0 || *n == 0) {
*scale = 1.;
return;
}
// Clean char * arguments
const char cleantranA = notransA ? 'N' : 'C';
const char cleantranB = notransB ? 'N' : 'C';

4
relapack/src/ztrtri.c
View File

@ -69,8 +69,8 @@ static void RELAPACK_ztrtri_rec(
}
// Constants
const double ONE[] = { 1. };
const double MONE[] = { -1. };
const double ONE[] = { 1., 0. };
const double MONE[] = { -1. , 0. };
// Splitting
const blasint n1 = ZREC_SPLIT(*n);