From d58cae90b98fda6d9e09d26e314f3db0bb1aa5da Mon Sep 17 00:00:00 2001
From: Daniel Waggoner <dwaggoner@frbatlanta.org>
Date: Thu, 3 Mar 2011 22:39:50 -0500
Subject: [PATCH] finished changes in bmatrix to make compatible with
blas/lapack calls in dynare (cherry picked from commit
390bb6746ea7c9a1f2eea11d3b420661789c9ebb)
---
include/blas_lapack.h | 13 ++++++-------
matrix/bmatrix.c | 39 +++++++++++++++++++++++++--------------
2 files changed, 31 insertions(+), 21 deletions(-)
diff --git a/include/blas_lapack.h b/include/blas_lapack.h
index 186cfc2..08ad286 100644
--- a/include/blas_lapack.h
+++ b/include/blas_lapack.h
@@ -18,7 +18,6 @@
#ifndef __BLAS_LAPACK__
#define __BLAS_LAPACK__
-
#if !defined(MATLAB_MEX_FILE) && !defined(OCTAVE_MEX_FILE)
#include "dw_std.h"
@@ -159,14 +158,14 @@ void sgelqf(lapack_int*,lapack_int*,float*,lapack_int*,float*,float*,lapack_int*
void dorglq(lapack_int*,lapack_int*,lapack_int*,double*,lapack_int*,double*,double*,lapack_int*,lapack_int*);
void sorglq(lapack_int*,lapack_int*,lapack_int*,float*,lapack_int*,float*,float*,lapack_int*,lapack_int*);
-void dgges(char*,char*,char*,void*,lapack_int*,double*,lapack_int*,double*,lapack_int*,lapack_int*,double*,double*,double*,double*,lapack_int*,double*,lapack_int*,double*,lapack_int*,void*,lapack_int*);
-void sgges(char*,char*,char*,void*,lapack_int*,float*,lapack_int*,float*,lapack_int*,lapack_int*,float*,float*,float*,float*,lapack_int*,float*,lapack_int*,float*,lapack_int*,void*,lapack_int*);
+void dgges(char*,char*,char*,lapack_int*,lapack_int*,double*,lapack_int*,double*,lapack_int*,lapack_int*,double*,double*,double*,double*,lapack_int*,double*,lapack_int*,double*,lapack_int*,void*,lapack_int*);
+void sgges(char*,char*,char*,lapack_int*,lapack_int*,float*,lapack_int*,float*,lapack_int*,lapack_int*,float*,float*,float*,float*,lapack_int*,float*,lapack_int*,float*,lapack_int*,void*,lapack_int*);
-void dtgsen(lapack_int*,void*,void*,void*,lapack_int*,double*,lapack_int*,double*,lapack_int*,double*,double*,double*,double*,lapack_int*,double*,lapack_int*,lapack_int*,double*,double*,double*,double*,lapack_int*,lapack_int*,lapack_int*,lapack_int*);
-void stgsen(lapack_int*,void*,void*,void*,lapack_int*,float*,lapack_int*,float*,lapack_int*,float*,float*,float*,float*,lapack_int*,float*,lapack_int*,lapack_int*,float*,float*,float*,float*,lapack_int*,lapack_int*,lapack_int*,lapack_int*);
+void dtgsen(lapack_int*,lapack_int*,lapack_int*,lapack_int*,lapack_int*,double*,lapack_int*,double*,lapack_int*,double*,double*,double*,double*,lapack_int*,double*,lapack_int*,lapack_int*,double*,double*,double*,double*,lapack_int*,lapack_int*,lapack_int*,lapack_int*);
+void stgsen(lapack_int*,lapack_int*,lapack_int*,lapack_int*,lapack_int*,float*,lapack_int*,float*,lapack_int*,float*,float*,float*,float*,lapack_int*,float*,lapack_int*,lapack_int*,float*,float*,float*,float*,lapack_int*,lapack_int*,lapack_int*,lapack_int*);
-void dtgexc(void*,void*,lapack_int*,double*,lapack_int*,double*,lapack_int*,double*,lapack_int*,double*,lapack_int*,lapack_int*,lapack_int*,double*,lapack_int*,lapack_int*);
-void stgexc(void*,void*,lapack_int*,float*,lapack_int*,float*,lapack_int*,float*,lapack_int*,float*,lapack_int*,lapack_int*,lapack_int*,float*,lapack_int*,lapack_int*);
+void dtgexc(lapack_int*,lapack_int*,lapack_int*,double*,lapack_int*,double*,lapack_int*,double*,lapack_int*,double*,lapack_int*,lapack_int*,lapack_int*,double*,lapack_int*,lapack_int*);
+void stgexc(lapack_int*,lapack_int*,lapack_int*,float*,lapack_int*,float*,lapack_int*,float*,lapack_int*,float*,lapack_int*,lapack_int*,lapack_int*,float*,lapack_int*,lapack_int*);
void dsyev(char*,char*,lapack_int*,double*,lapack_int*,double*,double*,lapack_int*,lapack_int*);
void ssyev(char*,char*,lapack_int*,float*,lapack_int*,float*,float*,lapack_int*,lapack_int*);
diff --git a/matrix/bmatrix.c b/matrix/bmatrix.c
index d7d5146..bc8e9e8 100644
--- a/matrix/bmatrix.c
+++ b/matrix/bmatrix.c
@@ -1465,13 +1465,13 @@ int bGeneralizedSchur_real(PRECISION *Q, PRECISION *Z, PRECISION *S, PRECISION *
if (!bt) bTransposeInPlace(B,n);
lwork=-1;
- gges(&jobvsl,&jobvsr,&sort,(void*)NULL,&bn,S,&bn,T,&bn,&simd,palpha_r,palpha_i,pbeta,Q,&bn,Z,&bn,&size,&lwork,(void*)NULL,&info);
+ gges(&jobvsl,&jobvsr,&sort,(lapack_int*)NULL,&bn,S,&bn,T,&bn,&simd,palpha_r,palpha_i,pbeta,Q,&bn,Z,&bn,&size,&lwork,(void*)NULL,&info);
if (!info)
if (!(work=(PRECISION*)dw_malloc((lwork=(lapack_int)size)*sizeof(PRECISION))))
rtrn=MEM_ERR;
else
{
- gges(&jobvsl,&jobvsr,&sort,(void*)NULL,&bn,S,&bn,T,&bn,&simd,palpha_r,palpha_i,pbeta,Q,&bn,Z,&bn,work,&lwork,(void*)NULL,&info);
+ gges(&jobvsl,&jobvsr,&sort,(lapack_int*)NULL,&bn,S,&bn,T,&bn,&simd,palpha_r,palpha_i,pbeta,Q,&bn,Z,&bn,work,&lwork,(void*)NULL,&info);
if (!info)
{
if (Q && !qt) bTransposeInPlace(Q,n);
@@ -1560,7 +1560,8 @@ int bReorderGeneralizedSchur_real(int *select, PRECISION *QQ, PRECISION *ZZ, PRE
#define tgsen dtgsen
#endif
- int ijob=0, wantq, wantz, lwork=-1, liwork=1, m, info, rtrn=NO_ERR, iwork;
+ int i, rtrn=NO_ERR;
+ lapack_int ijob=0, wantq, wantz, lwork=-1, liwork=-1, isize, info, *iwork, *bselect, bm, bn=n;
PRECISION size, *palpha_r, *palpha_i, *pbeta, *work, pl, pr, dif, *pQQ=QQ;
wantq=(QQ && Q) ? 1 : 0;
@@ -1569,8 +1570,9 @@ int bReorderGeneralizedSchur_real(int *select, PRECISION *QQ, PRECISION *ZZ, PRE
palpha_r=alpha_r ? alpha_r : (PRECISION*)dw_malloc(n*sizeof(PRECISION));
palpha_i=alpha_i ? alpha_i : (PRECISION*)dw_malloc(n*sizeof(PRECISION));
pbeta=beta ? beta : (PRECISION*)dw_malloc(n*sizeof(PRECISION));
+ bselect=(lapack_int*)dw_malloc(n*sizeof(lapack_int));
- if (palpha_r && palpha_i && pbeta)
+ if (palpha_r && palpha_i && pbeta && bselect)
{
if (SS != S)
if (st)
@@ -1610,18 +1612,25 @@ int bReorderGeneralizedSchur_real(int *select, PRECISION *QQ, PRECISION *ZZ, PRE
else
if (!zt) bTransposeInPlace(Z,n);
+ for (i=n-1; i >= 0; i--) bselect[i]=select[i];
- tgsen(&ijob,&wantq,&wantz,select,&n,SS,&n,TT,&n,palpha_r,palpha_i,pbeta,pQQ,&n,ZZ,&n,&m,
- &pl,&pr,&dif,&size,&lwork,&iwork,&liwork,&info);
+ tgsen(&ijob,&wantq,&wantz,bselect,&bn,SS,&bn,TT,&bn,palpha_r,palpha_i,pbeta,pQQ,&n,ZZ,&bn,&bm,
+ &pl,&pr,&dif,&size,&lwork,&isize,&liwork,&info);
if (!info)
if (!(work=(PRECISION*)dw_malloc((lwork=(int)size)*sizeof(PRECISION))))
rtrn=MEM_ERR;
else
{
- tgsen(&ijob,&wantq,&wantz,select,&n,SS,&n,TT,&n,palpha_r,palpha_i,pbeta,pQQ,&n,ZZ,&n,&m,
- &pl,&pr,&dif,work,&lwork,&iwork,&liwork,&info);
- if (info)
- rtrn=BLAS_LAPACK_ERR;
+ if (!(iwork=(lapack_int*)dw_malloc((liwork=isize)*sizeof(lapack_int))))
+ rtrn=MEM_ERR;
+ else
+ {
+ tgsen(&ijob,&wantq,&wantz,bselect,&bn,SS,&bn,TT,&bn,palpha_r,palpha_i,pbeta,pQQ,&bn,ZZ,&bn,&bm,
+ &pl,&pr,&dif,work,&lwork,iwork,&liwork,&info);
+ if (info)
+ rtrn=BLAS_LAPACK_ERR;
+ dw_free(iwork);
+ }
dw_free(work);
}
else
@@ -1640,6 +1649,7 @@ int bReorderGeneralizedSchur_real(int *select, PRECISION *QQ, PRECISION *ZZ, PRE
if (!alpha_r && palpha_r) dw_free(palpha_r);
if (!alpha_i && palpha_i) dw_free(palpha_i);
if (!beta && pbeta) dw_free(pbeta);
+ if (bselect) dw_free(bselect);
return rtrn;
@@ -1708,7 +1718,8 @@ int bSortGeneralizedSchur_real(PRECISION *QQ, PRECISION *ZZ, PRECISION *SS, PREC
#define tgexc dtgexc
#endif
- int wantq, wantz, lwork=-1, info, rtrn=NO_ERR, k, i, j, ii=2, jj=1;
+ int rtrn=NO_ERR, k, i, j;
+ lapack_int wantq, wantz, lwork=-1, info, ii=2, jj=1, bn=n;
PRECISION tmp1, tmp2, size, *work;
wantq=(QQ && Q) ? 1 : 0;
@@ -1750,9 +1761,9 @@ int bSortGeneralizedSchur_real(PRECISION *QQ, PRECISION *ZZ, PRECISION *SS, PREC
if (n > 1)
{
- tgexc(&wantq,&wantz,&n,SS,&n,TT,&n,QQ,&n,ZZ,&n,&ii,&jj,&size,&lwork,&info);
+ tgexc(&wantq,&wantz,&bn,SS,&bn,TT,&bn,QQ,&bn,ZZ,&bn,&ii,&jj,&size,&lwork,&info);
if (!info)
- if (!(work=(PRECISION*)dw_malloc((lwork=(int)size)*sizeof(PRECISION))))
+ if (!(work=(PRECISION*)dw_malloc((lwork=(lapack_int)size)*sizeof(PRECISION))))
rtrn=MEM_ERR;
else
{
@@ -1770,7 +1781,7 @@ int bSortGeneralizedSchur_real(PRECISION *QQ, PRECISION *ZZ, PRECISION *SS, PREC
k=((i == n-1) || (SS[i*n+i+1] == 0.0)) ? 1 : 2;
if (jj < ii)
{
- tgexc(&wantq,&wantz,&n,SS,&n,TT,&n,QQ,&n,ZZ,&n,&ii,&jj,work,&lwork,&info);
+ tgexc(&wantq,&wantz,&bn,SS,&bn,TT,&bn,QQ,&bn,ZZ,&bn,&ii,&jj,work,&lwork,&info);
if (info)
{
dw_free(work);
--
GitLab