kronecker DLLs: various modernizations and simplifications
parent
de159c0480
commit
1199d4abae
|
@ -1,5 +1,5 @@
|
|||
/*
|
||||
* Copyright © 2007-2011 Dynare Team
|
||||
* Copyright © 2007-2019 Dynare Team
|
||||
*
|
||||
* This file is part of Dynare.
|
||||
*
|
||||
|
@ -22,8 +22,6 @@
|
|||
* one can consider large matrices B and/or C.
|
||||
*/
|
||||
|
||||
#include <string.h>
|
||||
|
||||
#include <dynmex.h>
|
||||
#include <dynblas.h>
|
||||
|
||||
|
@ -34,10 +32,10 @@
|
|||
#define DEBUG_OMP 0
|
||||
|
||||
void
|
||||
full_A_times_kronecker_B_C(double *A, double *B, double *C, double *D,
|
||||
full_A_times_kronecker_B_C(const double *A, const double *B, const double *C, double *D,
|
||||
blas_int mA, blas_int nA, blas_int mB, blas_int nB, blas_int mC, blas_int nC, int number_of_threads)
|
||||
{
|
||||
#if USE_OMP
|
||||
#ifdef USE_OMP
|
||||
# pragma omp parallel for num_threads(number_of_threads)
|
||||
for (blas_int colD = 0; colD < nB*nC; colD++)
|
||||
{
|
||||
|
@ -54,23 +52,20 @@ full_A_times_kronecker_B_C(double *A, double *B, double *C, double *D,
|
|||
blas_int idxD = colD*mA;
|
||||
double BC = B[colB*mB+rowB]*C[colC*mC+rowC];
|
||||
for (blas_int rowD = 0; rowD < mA; rowD++)
|
||||
{
|
||||
D[idxD+rowD] += A[idxA+rowD]*BC;
|
||||
}
|
||||
D[idxD+rowD] += A[idxA+rowD]*BC;
|
||||
}
|
||||
}
|
||||
#else
|
||||
const blas_int shiftA = mA*mC;
|
||||
const blas_int shiftD = mA*nC;
|
||||
blas_int kd = 0, ka = 0;
|
||||
char transpose[2] = "N";
|
||||
double one = 1.0;
|
||||
for (blas_int col = 0; col < nB; col++)
|
||||
{
|
||||
ka = 0;
|
||||
for (blas_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);
|
||||
dgemm("N", "N", &mA, &nC, &mC, &B[mB*col+row], &A[ka], &mA, C, &mC, &one, &D[kd], &mA);
|
||||
ka += shiftA;
|
||||
}
|
||||
kd += shiftD;
|
||||
|
@ -79,9 +74,9 @@ full_A_times_kronecker_B_C(double *A, double *B, double *C, double *D,
|
|||
}
|
||||
|
||||
void
|
||||
full_A_times_kronecker_B_B(double *A, double *B, double *D, blas_int mA, blas_int nA, blas_int mB, blas_int nB, int number_of_threads)
|
||||
full_A_times_kronecker_B_B(const double *A, const double *B, double *D, blas_int mA, blas_int nA, blas_int mB, blas_int nB, int number_of_threads)
|
||||
{
|
||||
#if USE_OMP
|
||||
#ifdef USE_OMP
|
||||
# pragma omp parallel for num_threads(number_of_threads)
|
||||
for (blas_int colD = 0; colD < nB*nB; colD++)
|
||||
{
|
||||
|
@ -98,23 +93,20 @@ full_A_times_kronecker_B_B(double *A, double *B, double *D, blas_int mA, blas_in
|
|||
blas_int idxD = colD*mA;
|
||||
double BB = B[j1B*mB+i1B]*B[j2B*mB+i2B];
|
||||
for (blas_int rowD = 0; rowD < mA; rowD++)
|
||||
{
|
||||
D[idxD+rowD] += A[idxA+rowD]*BB;
|
||||
}
|
||||
D[idxD+rowD] += A[idxA+rowD]*BB;
|
||||
}
|
||||
}
|
||||
#else
|
||||
const blas_int shiftA = mA*mB;
|
||||
const blas_int shiftD = mA*nB;
|
||||
blas_int kd = 0, ka = 0;
|
||||
char transpose[2] = "N";
|
||||
double one = 1.0;
|
||||
for (blas_int col = 0; col < nB; col++)
|
||||
{
|
||||
ka = 0;
|
||||
for (blas_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);
|
||||
dgemm("N", "N", &mA, &nB, &mB, &B[mB*col+row], &A[ka], &mA, B, &mB, &one, &D[kd], &mA);
|
||||
ka += shiftA;
|
||||
}
|
||||
kd += shiftD;
|
||||
|
@ -130,11 +122,11 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
|||
DYN_MEX_FUNC_ERR_MSG_TXT("A_times_B_kronecker_C takes 3 or 4 input arguments and provides 2 output arguments.");
|
||||
|
||||
// Get & Check dimensions (columns and rows):
|
||||
mwSize mA, nA, mB, nB, mC, nC;
|
||||
mA = mxGetM(prhs[0]);
|
||||
nA = mxGetN(prhs[0]);
|
||||
mB = mxGetM(prhs[1]);
|
||||
nB = mxGetN(prhs[1]);
|
||||
size_t mA = mxGetM(prhs[0]);
|
||||
size_t nA = mxGetN(prhs[0]);
|
||||
size_t mB = mxGetM(prhs[1]);
|
||||
size_t nB = mxGetN(prhs[1]);
|
||||
size_t mC, nC;
|
||||
if (nrhs == 4) // A*kron(B,C) is to be computed.
|
||||
{
|
||||
mC = mxGetM(prhs[2]);
|
||||
|
@ -148,10 +140,10 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
|||
DYN_MEX_FUNC_ERR_MSG_TXT("Input dimension error!");
|
||||
}
|
||||
// Get input matrices:
|
||||
double *B, *C, *A;
|
||||
int numthreads;
|
||||
A = mxGetPr(prhs[0]);
|
||||
B = mxGetPr(prhs[1]);
|
||||
const double *A = mxGetPr(prhs[0]);
|
||||
const double *B = mxGetPr(prhs[1]);
|
||||
const double *C;
|
||||
if (nrhs == 4)
|
||||
{
|
||||
C = mxGetPr(prhs[2]);
|
||||
|
@ -161,24 +153,17 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
|||
numthreads = static_cast<int>(mxGetScalar(prhs[2]));
|
||||
|
||||
// Initialization of the ouput:
|
||||
double *D;
|
||||
if (nrhs == 4)
|
||||
{
|
||||
plhs[0] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
|
||||
}
|
||||
plhs[0] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
|
||||
else
|
||||
{
|
||||
plhs[0] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
|
||||
}
|
||||
D = mxGetPr(plhs[0]);
|
||||
plhs[0] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
|
||||
double *D = mxGetPr(plhs[0]);
|
||||
|
||||
// Computational part:
|
||||
if (nrhs == 3)
|
||||
{
|
||||
full_A_times_kronecker_B_B(A, B, &D[0], mA, nA, mB, nB, numthreads);
|
||||
}
|
||||
full_A_times_kronecker_B_B(A, B, D, mA, nA, mB, nB, numthreads);
|
||||
else
|
||||
{
|
||||
full_A_times_kronecker_B_C(A, B, C, &D[0], mA, nA, mB, nB, mC, nC, numthreads);
|
||||
}
|
||||
full_A_times_kronecker_B_C(A, B, C, D, mA, nA, mB, nB, mC, nC, numthreads);
|
||||
|
||||
plhs[1] = mxCreateDoubleScalar(0);
|
||||
}
|
||||
|
|
|
@ -1,5 +1,5 @@
|
|||
/*
|
||||
* Copyright © 2007-2017 Dynare Team
|
||||
* Copyright © 2007-2019 Dynare Team
|
||||
*
|
||||
* This file is part of Dynare.
|
||||
*
|
||||
|
@ -23,7 +23,7 @@
|
|||
* (dynare format). This mex file should not be used outside dr1.m.
|
||||
*/
|
||||
|
||||
#include <string.h>
|
||||
#include <algorithm>
|
||||
|
||||
#include <dynmex.h>
|
||||
|
||||
|
@ -34,8 +34,8 @@
|
|||
#define DEBUG_OMP 0
|
||||
|
||||
void
|
||||
sparse_hessian_times_B_kronecker_B(const mwIndex *isparseA, const mwIndex *jsparseA, double *vsparseA,
|
||||
double *B, double *D, mwSize mA, mwSize nA, mwSize mB, mwSize nB, int number_of_threads)
|
||||
sparse_hessian_times_B_kronecker_B(const mwIndex *isparseA, const mwIndex *jsparseA, const double *vsparseA,
|
||||
const double *B, double *D, size_t mA, size_t nA, size_t mB, size_t nB, int number_of_threads)
|
||||
{
|
||||
/*
|
||||
** Loop over the columns of kron(B,B) (or of the result matrix D).
|
||||
|
@ -45,12 +45,12 @@ sparse_hessian_times_B_kronecker_B(const mwIndex *isparseA, const mwIndex *jspar
|
|||
#if USE_OMP
|
||||
# pragma omp parallel for num_threads(number_of_threads)
|
||||
#endif
|
||||
for (mwIndex j1B = 0; j1B < nB; j1B++)
|
||||
for (mwIndex j1B = 0; j1B < static_cast<mwIndex>(nB); j1B++)
|
||||
{
|
||||
#if DEBUG_OMP
|
||||
mexPrintf("%d thread number is %d (%d).\n", j1B, omp_get_thread_num(), omp_get_num_threads());
|
||||
#endif
|
||||
for (mwIndex j2B = j1B; j2B < nB; j2B++)
|
||||
for (mwIndex j2B = j1B; j2B < static_cast<mwIndex>(nB); j2B++)
|
||||
{
|
||||
mwIndex jj = j1B*nB+j2B; // column of kron(B,B) index.
|
||||
mwIndex iv = 0;
|
||||
|
@ -60,16 +60,16 @@ sparse_hessian_times_B_kronecker_B(const mwIndex *isparseA, const mwIndex *jspar
|
|||
/*
|
||||
** Loop over the rows of kron(B,B) (column jj).
|
||||
*/
|
||||
for (mwIndex ii = 0; ii < nA; ii++)
|
||||
for (mwIndex ii = 0; ii < static_cast<mwIndex>(nA); ii++)
|
||||
{
|
||||
k1 = jsparseA[ii];
|
||||
k2 = jsparseA[ii+1];
|
||||
if (k1 < k2) // otherwise column ii of A does not have non zero elements (and there is nothing to compute).
|
||||
{
|
||||
++nz_in_column_ii_of_A;
|
||||
mwIndex i1B = (ii/mB);
|
||||
mwIndex i2B = (ii%mB);
|
||||
double bb = B[j1B*mB+i1B]*B[j2B*mB+i2B];
|
||||
mwIndex i1B = ii / mB;
|
||||
mwIndex i2B = ii % mB;
|
||||
double bb = B[j1B*mB+i1B]*B[j2B*mB+i2B];
|
||||
/*
|
||||
** Loop over the non zero entries of A(:,ii).
|
||||
*/
|
||||
|
@ -82,17 +82,15 @@ sparse_hessian_times_B_kronecker_B(const mwIndex *isparseA, const mwIndex *jspar
|
|||
}
|
||||
}
|
||||
if (nz_in_column_ii_of_A > 0)
|
||||
{
|
||||
memcpy(&D[(j2B*nB+j1B)*mA], &D[jj*mA], mA*sizeof(double));
|
||||
}
|
||||
std::copy_n(&D[jj*mA], mA, &D[(j2B*nB+j1B)*mA]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
sparse_hessian_times_B_kronecker_C(const mwIndex *isparseA, const mwIndex *jsparseA, double *vsparseA,
|
||||
double *B, double *C, double *D,
|
||||
mwSize mA, mwSize nA, mwSize mB, mwSize nB, mwSize mC, mwSize nC, int number_of_threads)
|
||||
sparse_hessian_times_B_kronecker_C(const mwIndex *isparseA, const mwIndex *jsparseA, const double *vsparseA,
|
||||
const double *B, const double *C, double *D,
|
||||
size_t mA, size_t nA, size_t mB, size_t nB, size_t mC, size_t nC, int number_of_threads)
|
||||
{
|
||||
/*
|
||||
** Loop over the columns of kron(B,B) (or of the result matrix D).
|
||||
|
@ -100,7 +98,7 @@ sparse_hessian_times_B_kronecker_C(const mwIndex *isparseA, const mwIndex *jspar
|
|||
#if USE_OMP
|
||||
# pragma omp parallel for num_threads(number_of_threads)
|
||||
#endif
|
||||
for (mwIndex jj = 0; jj < nB*nC; jj++) // column of kron(B,C) index.
|
||||
for (mwIndex jj = 0; jj < static_cast<mwIndex>(nB*nC); jj++) // column of kron(B,C) index.
|
||||
{
|
||||
// Uncomment the following line to check if all processors are used.
|
||||
#if DEBUG_OMP
|
||||
|
@ -115,15 +113,15 @@ sparse_hessian_times_B_kronecker_C(const mwIndex *isparseA, const mwIndex *jspar
|
|||
/*
|
||||
** Loop over the rows of kron(B,C) (column jj).
|
||||
*/
|
||||
for (mwIndex ii = 0; ii < nA; ii++)
|
||||
for (mwIndex ii = 0; ii < static_cast<mwIndex>(nA); ii++)
|
||||
{
|
||||
k1 = jsparseA[ii];
|
||||
k2 = jsparseA[ii+1];
|
||||
if (k1 < k2) // otherwise column ii of A does not have non zero elements (and there is nothing to compute).
|
||||
{
|
||||
++nz_in_column_ii_of_A;
|
||||
mwIndex iC = (ii%mB);
|
||||
mwIndex iB = (ii/mB);
|
||||
mwIndex iC = ii % mB;
|
||||
mwIndex iB = ii / mB;
|
||||
double cb = C[jC*mC+iC]*B[jB*mB+iB];
|
||||
/*
|
||||
** Loop over the non zero entries of A(:,ii).
|
||||
|
@ -143,18 +141,18 @@ void
|
|||
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
||||
{
|
||||
// Check input and output:
|
||||
if ((nrhs > 4) || (nrhs < 3))
|
||||
if (nrhs > 4 || nrhs < 3)
|
||||
DYN_MEX_FUNC_ERR_MSG_TXT("sparse_hessian_times_B_kronecker_C takes 3 or 4 input arguments and provides 2 output arguments.");
|
||||
|
||||
if (!mxIsSparse(prhs[0]))
|
||||
DYN_MEX_FUNC_ERR_MSG_TXT("sparse_hessian_times_B_kronecker_C: First input must be a sparse (dynare) hessian matrix.");
|
||||
|
||||
// Get & Check dimensions (columns and rows):
|
||||
mwSize mA, nA, mB, nB, mC, nC;
|
||||
mA = mxGetM(prhs[0]);
|
||||
nA = mxGetN(prhs[0]);
|
||||
mB = mxGetM(prhs[1]);
|
||||
nB = mxGetN(prhs[1]);
|
||||
size_t mA = mxGetM(prhs[0]);
|
||||
size_t nA = mxGetN(prhs[0]);
|
||||
size_t mB = mxGetM(prhs[1]);
|
||||
size_t nB = mxGetN(prhs[1]);
|
||||
size_t mC, nC;
|
||||
if (nrhs == 4) // A*kron(B,C) is to be computed.
|
||||
{
|
||||
mC = mxGetM(prhs[2]);
|
||||
|
@ -168,9 +166,9 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
|||
DYN_MEX_FUNC_ERR_MSG_TXT("Input dimension error!");
|
||||
}
|
||||
// Get input matrices:
|
||||
double *B, *C;
|
||||
int numthreads;
|
||||
B = mxGetPr(prhs[1]);
|
||||
const double *B = mxGetPr(prhs[1]);
|
||||
const double *C;
|
||||
numthreads = static_cast<int>(mxGetScalar(prhs[2]));
|
||||
if (nrhs == 4)
|
||||
{
|
||||
|
@ -180,26 +178,19 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
|||
// Sparse (dynare) hessian matrix.
|
||||
const mwIndex *isparseA = mxGetIr(prhs[0]);
|
||||
const mwIndex *jsparseA = mxGetJc(prhs[0]);
|
||||
double *vsparseA = mxGetPr(prhs[0]);
|
||||
const double *vsparseA = mxGetPr(prhs[0]);
|
||||
// Initialization of the ouput:
|
||||
double *D;
|
||||
if (nrhs == 4)
|
||||
{
|
||||
plhs[0] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
|
||||
}
|
||||
plhs[0] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
|
||||
else
|
||||
{
|
||||
plhs[0] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
|
||||
}
|
||||
D = mxGetPr(plhs[0]);
|
||||
plhs[0] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
|
||||
double *D = mxGetPr(plhs[0]);
|
||||
|
||||
// Computational part:
|
||||
if (nrhs == 3)
|
||||
{
|
||||
sparse_hessian_times_B_kronecker_B(isparseA, jsparseA, vsparseA, B, D, mA, nA, mB, nB, numthreads);
|
||||
}
|
||||
sparse_hessian_times_B_kronecker_B(isparseA, jsparseA, vsparseA, B, D, mA, nA, mB, nB, numthreads);
|
||||
else
|
||||
{
|
||||
sparse_hessian_times_B_kronecker_C(isparseA, jsparseA, vsparseA, B, C, D, mA, nA, mB, nB, mC, nC, numthreads);
|
||||
}
|
||||
sparse_hessian_times_B_kronecker_C(isparseA, jsparseA, vsparseA, B, C, D, mA, nA, mB, nB, mC, nC, numthreads);
|
||||
|
||||
plhs[1] = mxCreateDoubleScalar(0);
|
||||
}
|
||||
|
|
Loading…
Reference in New Issue