Fix integer types in Kronecker DLLs
parent
23d797d214
commit
ffc53ca910
|
@ -1,5 +1,5 @@
|
||||||
/*
|
/*
|
||||||
* Copyright (C) 2007-2010 Dynare Team
|
* Copyright (C) 2007-2011 Dynare Team
|
||||||
*
|
*
|
||||||
* This file is part of Dynare.
|
* This file is part of Dynare.
|
||||||
*
|
*
|
||||||
|
@ -39,36 +39,36 @@ full_A_times_kronecker_B_C(double *A, double *B, double *C, double *D,
|
||||||
{
|
{
|
||||||
#if USE_OMP
|
#if USE_OMP
|
||||||
# pragma omp parallel for num_threads(number_of_threads)
|
# pragma omp parallel for num_threads(number_of_threads)
|
||||||
for (long int colD = 0; colD < nB*nC; colD++)
|
for (blas_int colD = 0; colD < nB*nC; colD++)
|
||||||
{
|
{
|
||||||
# if DEBUG_OMP
|
# if DEBUG_OMP
|
||||||
mexPrintf("%d thread number is %d (%d).\n", colD, omp_get_thread_num(), omp_get_num_threads());
|
mexPrintf("%d thread number is %d (%d).\n", colD, omp_get_thread_num(), omp_get_num_threads());
|
||||||
# endif
|
# endif
|
||||||
unsigned int colB = colD/nC;
|
blas_int colB = colD/nC;
|
||||||
unsigned int colC = colD%nC;
|
blas_int colC = colD%nC;
|
||||||
for (int colA = 0; colA < nA; colA++)
|
for (blas_int colA = 0; colA < nA; colA++)
|
||||||
{
|
{
|
||||||
unsigned int rowB = colA/mC;
|
blas_int rowB = colA/mC;
|
||||||
unsigned int rowC = colA%mC;
|
blas_int rowC = colA%mC;
|
||||||
unsigned int idxA = colA*mA;
|
blas_int idxA = colA*mA;
|
||||||
unsigned int idxD = colD*mA;
|
blas_int idxD = colD*mA;
|
||||||
double BC = B[colB*mB+rowB]*C[colC*mC+rowC];
|
double BC = B[colB*mB+rowB]*C[colC*mC+rowC];
|
||||||
for (int rowD = 0; rowD < mA; rowD++)
|
for (blas_int rowD = 0; rowD < mA; rowD++)
|
||||||
{
|
{
|
||||||
D[idxD+rowD] += A[idxA+rowD]*BC;
|
D[idxD+rowD] += A[idxA+rowD]*BC;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
#else
|
#else
|
||||||
const unsigned long shiftA = mA*mC;
|
const blas_int shiftA = mA*mC;
|
||||||
const unsigned long shiftD = mA*nC;
|
const blas_int shiftD = mA*nC;
|
||||||
unsigned long int kd = 0, ka = 0;
|
blas_int kd = 0, ka = 0;
|
||||||
char transpose[2] = "N";
|
char transpose[2] = "N";
|
||||||
double one = 1.0;
|
double one = 1.0;
|
||||||
for (unsigned long int col = 0; col < nB; col++)
|
for (blas_int col = 0; col < nB; col++)
|
||||||
{
|
{
|
||||||
ka = 0;
|
ka = 0;
|
||||||
for (unsigned long int row = 0; row < mB; row++)
|
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(transpose, transpose, &mA, &nC, &mC, &B[mB*col+row], &A[ka], &mA, &C[0], &mC, &one, &D[kd], &mA);
|
||||||
ka += shiftA;
|
ka += shiftA;
|
||||||
|
@ -83,36 +83,36 @@ full_A_times_kronecker_B_B(double *A, double *B, double *D, blas_int mA, blas_in
|
||||||
{
|
{
|
||||||
#if USE_OMP
|
#if USE_OMP
|
||||||
# pragma omp parallel for num_threads(number_of_threads)
|
# pragma omp parallel for num_threads(number_of_threads)
|
||||||
for (long int colD = 0; colD < nB*nB; colD++)
|
for (blas_int colD = 0; colD < nB*nB; colD++)
|
||||||
{
|
{
|
||||||
# if DEBUG_OMP
|
# if DEBUG_OMP
|
||||||
mexPrintf("%d thread number is %d (%d).\n", colD, omp_get_thread_num(), omp_get_num_threads());
|
mexPrintf("%d thread number is %d (%d).\n", colD, omp_get_thread_num(), omp_get_num_threads());
|
||||||
# endif
|
# endif
|
||||||
unsigned int j1B = colD/nB;
|
blas_int j1B = colD/nB;
|
||||||
unsigned int j2B = colD%nB;
|
blas_int j2B = colD%nB;
|
||||||
for (int colA = 0; colA < nA; colA++)
|
for (blas_int colA = 0; colA < nA; colA++)
|
||||||
{
|
{
|
||||||
unsigned int i1B = colA/mB;
|
blas_int i1B = colA/mB;
|
||||||
unsigned int i2B = colA%mB;
|
blas_int i2B = colA%mB;
|
||||||
unsigned int idxA = colA*mA;
|
blas_int idxA = colA*mA;
|
||||||
unsigned int idxD = colD*mA;
|
blas_int idxD = colD*mA;
|
||||||
double BB = B[j1B*mB+i1B]*B[j2B*mB+i2B];
|
double BB = B[j1B*mB+i1B]*B[j2B*mB+i2B];
|
||||||
for (int rowD = 0; rowD < mA; rowD++)
|
for (blas_int rowD = 0; rowD < mA; rowD++)
|
||||||
{
|
{
|
||||||
D[idxD+rowD] += A[idxA+rowD]*BB;
|
D[idxD+rowD] += A[idxA+rowD]*BB;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
#else
|
#else
|
||||||
const unsigned long int shiftA = mA*mB;
|
const blas_int shiftA = mA*mB;
|
||||||
const unsigned long int shiftD = mA*nB;
|
const blas_int shiftD = mA*nB;
|
||||||
unsigned long int kd = 0, ka = 0;
|
blas_int kd = 0, ka = 0;
|
||||||
char transpose[2] = "N";
|
char transpose[2] = "N";
|
||||||
double one = 1.0;
|
double one = 1.0;
|
||||||
for (unsigned long int col = 0; col < nB; col++)
|
for (blas_int col = 0; col < nB; col++)
|
||||||
{
|
{
|
||||||
ka = 0;
|
ka = 0;
|
||||||
for (unsigned long int row = 0; row < mB; row++)
|
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(transpose, transpose, &mA, &nB, &mB, &B[mB*col+row], &A[ka], &mA, &B[0], &mB, &one, &D[kd], &mA);
|
||||||
ka += shiftA;
|
ka += shiftA;
|
||||||
|
@ -174,11 +174,11 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
||||||
// Computational part:
|
// Computational part:
|
||||||
if (nrhs == 3)
|
if (nrhs == 3)
|
||||||
{
|
{
|
||||||
full_A_times_kronecker_B_B(A, B, &D[0], (int) mA, (int) nA, (int) mB, (int) nB, numthreads);
|
full_A_times_kronecker_B_B(A, B, &D[0], mA, nA, mB, nB, numthreads);
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
full_A_times_kronecker_B_C(A, B, C, &D[0], (int) mA, (int) nA, (int) mB, (int) nB, (int) mC, (int) nC, numthreads);
|
full_A_times_kronecker_B_C(A, B, C, &D[0], mA, nA, mB, nB, mC, nC, numthreads);
|
||||||
}
|
}
|
||||||
plhs[0] = mxCreateDoubleScalar(0);
|
plhs[0] = mxCreateDoubleScalar(0);
|
||||||
}
|
}
|
||||||
|
|
|
@ -1,5 +1,5 @@
|
||||||
/*
|
/*
|
||||||
* Copyright (C) 2007-2010 Dynare Team
|
* Copyright (C) 2007-2011 Dynare Team
|
||||||
*
|
*
|
||||||
* This file is part of Dynare.
|
* This file is part of Dynare.
|
||||||
*
|
*
|
||||||
|
@ -45,37 +45,37 @@ sparse_hessian_times_B_kronecker_B(mwIndex *isparseA, mwIndex *jsparseA, double
|
||||||
#if USE_OMP
|
#if USE_OMP
|
||||||
# pragma omp parallel for num_threads(number_of_threads)
|
# pragma omp parallel for num_threads(number_of_threads)
|
||||||
#endif
|
#endif
|
||||||
for (int j1B = 0; j1B < nB; j1B++)
|
for (mwIndex j1B = 0; j1B < nB; j1B++)
|
||||||
{
|
{
|
||||||
#if DEBUG_OMP
|
#if DEBUG_OMP
|
||||||
mexPrintf("%d thread number is %d (%d).\n", j1B, omp_get_thread_num(), omp_get_num_threads());
|
mexPrintf("%d thread number is %d (%d).\n", j1B, omp_get_thread_num(), omp_get_num_threads());
|
||||||
#endif
|
#endif
|
||||||
for (unsigned int j2B = j1B; j2B < nB; j2B++)
|
for (mwIndex j2B = j1B; j2B < nB; j2B++)
|
||||||
{
|
{
|
||||||
unsigned long int jj = j1B*nB+j2B; // column of kron(B,B) index.
|
mwIndex jj = j1B*nB+j2B; // column of kron(B,B) index.
|
||||||
unsigned long int iv = 0;
|
mwIndex iv = 0;
|
||||||
unsigned int nz_in_column_ii_of_A = 0;
|
int nz_in_column_ii_of_A = 0;
|
||||||
unsigned int k1 = 0;
|
mwIndex k1 = 0;
|
||||||
unsigned int k2 = 0;
|
mwIndex k2 = 0;
|
||||||
/*
|
/*
|
||||||
** Loop over the rows of kron(B,B) (column jj).
|
** Loop over the rows of kron(B,B) (column jj).
|
||||||
*/
|
*/
|
||||||
for (unsigned long int ii = 0; ii < nA; ii++)
|
for (mwIndex ii = 0; ii < nA; ii++)
|
||||||
{
|
{
|
||||||
k1 = jsparseA[ii];
|
k1 = jsparseA[ii];
|
||||||
k2 = jsparseA[ii+1];
|
k2 = jsparseA[ii+1];
|
||||||
if (k1 < k2) // otherwise column ii of A does not have non zero elements (and there is nothing to compute).
|
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;
|
++nz_in_column_ii_of_A;
|
||||||
unsigned int i1B = (ii/mB);
|
mwIndex i1B = (ii/mB);
|
||||||
unsigned int i2B = (ii%mB);
|
mwIndex i2B = (ii%mB);
|
||||||
double bb = B[j1B*mB+i1B]*B[j2B*mB+i2B];
|
double bb = B[j1B*mB+i1B]*B[j2B*mB+i2B];
|
||||||
/*
|
/*
|
||||||
** Loop over the non zero entries of A(:,ii).
|
** Loop over the non zero entries of A(:,ii).
|
||||||
*/
|
*/
|
||||||
for (unsigned int k = k1; k < k2; k++)
|
for (mwIndex k = k1; k < k2; k++)
|
||||||
{
|
{
|
||||||
unsigned int kk = isparseA[k];
|
mwIndex kk = isparseA[k];
|
||||||
D[jj*mA+kk] = D[jj*mA+kk] + bb*vsparseA[iv];
|
D[jj*mA+kk] = D[jj*mA+kk] + bb*vsparseA[iv];
|
||||||
iv++;
|
iv++;
|
||||||
}
|
}
|
||||||
|
@ -100,37 +100,37 @@ sparse_hessian_times_B_kronecker_C(mwIndex *isparseA, mwIndex *jsparseA, double
|
||||||
#if USE_OMP
|
#if USE_OMP
|
||||||
# pragma omp parallel for num_threads(number_of_threads)
|
# pragma omp parallel for num_threads(number_of_threads)
|
||||||
#endif
|
#endif
|
||||||
for (long int jj = 0; jj < nB*nC; jj++) // column of kron(B,C) index.
|
for (mwIndex jj = 0; jj < nB*nC; jj++) // column of kron(B,C) index.
|
||||||
{
|
{
|
||||||
// Uncomment the following line to check if all processors are used.
|
// Uncomment the following line to check if all processors are used.
|
||||||
#if DEBUG_OMP
|
#if DEBUG_OMP
|
||||||
mexPrintf("%d thread number is %d (%d).\n", jj, omp_get_thread_num(), omp_get_num_threads());
|
mexPrintf("%d thread number is %d (%d).\n", jj, omp_get_thread_num(), omp_get_num_threads());
|
||||||
#endif
|
#endif
|
||||||
unsigned int jB = jj/nC;
|
mwIndex jB = jj/nC;
|
||||||
unsigned int jC = jj%nC;
|
mwIndex jC = jj%nC;
|
||||||
unsigned int k1 = 0;
|
mwIndex k1 = 0;
|
||||||
unsigned int k2 = 0;
|
mwIndex k2 = 0;
|
||||||
unsigned long int iv = 0;
|
mwIndex iv = 0;
|
||||||
unsigned int nz_in_column_ii_of_A = 0;
|
int nz_in_column_ii_of_A = 0;
|
||||||
/*
|
/*
|
||||||
** Loop over the rows of kron(B,C) (column jj).
|
** Loop over the rows of kron(B,C) (column jj).
|
||||||
*/
|
*/
|
||||||
for (unsigned long int ii = 0; ii < nA; ii++)
|
for (mwIndex ii = 0; ii < nA; ii++)
|
||||||
{
|
{
|
||||||
k1 = jsparseA[ii];
|
k1 = jsparseA[ii];
|
||||||
k2 = jsparseA[ii+1];
|
k2 = jsparseA[ii+1];
|
||||||
if (k1 < k2) // otherwise column ii of A does not have non zero elements (and there is nothing to compute).
|
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;
|
++nz_in_column_ii_of_A;
|
||||||
unsigned int iC = (ii%mB);
|
mwIndex iC = (ii%mB);
|
||||||
unsigned int iB = (ii/mB);
|
mwIndex iB = (ii/mB);
|
||||||
double cb = C[jC*mC+iC]*B[jB*mB+iB];
|
double cb = C[jC*mC+iC]*B[jB*mB+iB];
|
||||||
/*
|
/*
|
||||||
** Loop over the non zero entries of A(:,ii).
|
** Loop over the non zero entries of A(:,ii).
|
||||||
*/
|
*/
|
||||||
for (unsigned int k = k1; k < k2; k++)
|
for (mwIndex k = k1; k < k2; k++)
|
||||||
{
|
{
|
||||||
unsigned int kk = isparseA[k];
|
mwIndex kk = isparseA[k];
|
||||||
D[jj*mA+kk] = D[jj*mA+kk] + cb*vsparseA[iv];
|
D[jj*mA+kk] = D[jj*mA+kk] + cb*vsparseA[iv];
|
||||||
iv++;
|
iv++;
|
||||||
}
|
}
|
||||||
|
|
Loading…
Reference in New Issue