diff --git a/mex/sources/kronecker/A_times_B_kronecker_C.cc b/mex/sources/kronecker/A_times_B_kronecker_C.cc index 5dc437ed9..72d631328 100644 --- a/mex/sources/kronecker/A_times_B_kronecker_C.cc +++ b/mex/sources/kronecker/A_times_B_kronecker_C.cc @@ -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 - #include #include @@ -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(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); } diff --git a/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc b/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc index 14e798032..1384f1dd8 100644 --- a/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc +++ b/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc @@ -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 +#include #include @@ -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(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(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(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(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(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(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); }