stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added mex file to compute A*kron(B,B) or A*kron(B,C). 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1443 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
adjemian 2007-11-15 14:58:56 +00:00
parent 8204b9abd6
commit d9db719b7b
2 changed files with 115 additions and 60 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

115
matlab/AkronBC.cc Normal file
View File

@ -0,0 +1,115 @@
// This mex file computes A*kron(B,C) or A*kron(B,B) without explicitely building kron(B,C) or kron(B,B), so that
// one can consider large matrices B and/or C.
//
// (linux)SYNTAX:
// mex AkronBC.cc /opt/matlab2007b/bin/glnx86/mkl.so
//
// stephane.adjemian@ens.fr [15-11-2007]
// Dynare Team, 2007.
#include "mex.h"
#include "blas.h"
void A_times_kronecker_B_B(double *A, double *B, double *D, int mA, int nA, int mB, int nB)
{
const int shiftA = mA*mB ;
const int shiftD = mA*nB ;
char transpose[2] = "N";
double one = 1.0;
int kd = 0 ;
for(int col=0; col<nB; col++)
{
int ka = 0 ;
for(int row=0; row<mB; row++)
{
dgemm(transpose, transpose, &mA, &nB, &mB, &B[mB*col+row], &A[ka], &mA, &B[0], &mB, &one, &D[kd], &mA);
ka += shiftA;
}
kd += shiftD;
}
}
void A_times_kronecker_B_C(double *A, double *B, double *C, double *D,
int mA, int nA, int mB, int nB, int mC, int nC)
{
const int shiftA = mA*mC ;
const int shiftD = mA*nC ;
double one = 1.0 ;
char transpose[2] = "N";
int kd = 0 ;
for(int col=0; col<nB; col++)
{
int ka = 0 ;
for(int row=0; row<mB; row++)
{
dgemm(transpose, transpose, &mA, &nC, &mC, &B[mB*col+row], &A[ka], &mA, &C[0], &mC, &one, &D[kd], &mA);
ka += shiftA;
}
kd += shiftD;
}
}
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{
// Check input and output:
if ( (nrhs > 3) || (nrhs <2) )
{
mexErrMsgTxt("Two or Three input arguments required.");
}
if (nlhs>1)
{
mexErrMsgTxt("Too many output arguments.");
}
// Get & Check dimensions (columns and rows):
int mA, nA, mB, nB, mC, nC;
mA = (int)mxGetM(prhs[0]);
nA = (int)mxGetN(prhs[0]);
mB = (int)mxGetM(prhs[1]);
nB = (int)mxGetN(prhs[1]);
if (nrhs == 3)// A*kron(B,C) is to be computed.
{
mC = (int)mxGetM(prhs[2]);
nC = (int)mxGetN(prhs[2]);
if (mB*mC != nA)
{
mexErrMsgTxt("Input dimension error!");
}
}
else// A*kron(B,B) is to be computed.
{
if (mB*mB != nA)
{
mexErrMsgTxt("Input dimension error!");
}
}
// Get input matrices:
double *A, *B, *C;
A = mxGetPr(prhs[0]);
B = mxGetPr(prhs[1]);
if (nrhs == 3)
{
C = mxGetPr(prhs[2]);
}
// Initialization of the ouput:
double *D;
if (nrhs == 3)
{
plhs[0] = mxCreateDoubleMatrix(mA,nB*nC,mxREAL);
}
else
{
plhs[0] = mxCreateDoubleMatrix(mA,nB*nB,mxREAL);
}
D = mxGetPr(plhs[0]);
// Computational part:
if (nrhs == 2)
{
A_times_kronecker_B_B(A, B, &D[0], mA, nA, mB, nB);
}
else
{
A_times_kronecker_B_C(A, B, C, &D[0], mA, nA, mB, nB, mC, nC);
}
}

60
matlab/prodk.cc
View File

@ -1,60 +0,0 @@
#include "mex.h"
#include "blas.h"
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{
int inc = 1;
char tA[2] = {'N','\0'};
double beta = 1.0;
// Check input and output:
if (nrhs != 3)
{
mexErrMsgTxt("Three input arguments required.");
}
else if (nlhs>1)
{
mexErrMsgTxt("Too many output arguments.");
}
// Get & Check dimensions:
int mA, nA, mb, nb, mc, nc;
mA = (int)mxGetM(prhs[0]);
nA = (int)mxGetN(prhs[0]);
mb = (int)mxGetM(prhs[1]);
nb = (int)mxGetN(prhs[1]);
if (nb != 1)
{
mexErrMsgTxt("Second argument must be a column vector.");
}
mc = (int)mxGetM(prhs[2]);
nc = (int)mxGetN(prhs[1]);
if (nc != 1)
{
mexErrMsgTxt("Third argument must be a column vector.");
}
if (mb*mc != nA)
{
mexErrMsgTxt("Input dimension error.");
}
// Get input matrices:
double *A, *b, *c ;
A = mxGetPr(prhs[0]);
b = mxGetPr(prhs[1]);
c = mxGetPr(prhs[2]);
// Initialization of the ouput:
double *d;
plhs[0] = mxCreateDoubleMatrix(mA,1,mxREAL);
d = mxGetPr(plhs[0]);
// Computational part:
int k = 0;
for (int i=0; i<mb ; i++)
{
dgemv(tA, &mA, &mc, &b[i], &A[k], &mA, &c[0], &inc, &beta, &d[0], &inc);
k += mc*mA;
}
}