/* The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2013 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU General Public License as published by the 
%   Free Software Foundation, either version 3 of the License, or (at your 
%   option) any later version.
%
% This program is distributed in the hope that it will be useful, but 
%   WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
%   General Public License for more details.
%
% You should have received a copy of the GNU General Public License along 
%   with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    /* Syntax: measureVector = evaluateMeasurePolynome(measureTensor, ...
     *           vecMatrix)
     */
    double *outMeas, *inTen, *inVec;
    mwSize nVec, nC, nDim;
    mwSize iVec, i, j, k;

    /* get sizes */
    nC = *(mxGetDimensions(prhs[1]));
    nVec = *(mxGetDimensions(prhs[1]) + 1);
    nDim = (mwSize) *mxGetPr(prhs[2]);

    /* create output vector*/
    plhs[0] = mxCreateDoubleMatrix(1, nVec, mxREAL);
    
    /* get data pointers */
    outMeas = mxGetPr(plhs[0]);
    inTen = mxGetPr(prhs[0]);
    inVec = mxGetPr(prhs[1]);
    
    /* cycle through vectors */
    for(iVec = 0; iVec < nVec; iVec++) {
        switch(nDim) {
            case 1:
                /* one-fold-loop */
                for(i = 0; i < nC; i++) {
                    *(outMeas + iVec) = *(outMeas + iVec) +
                            *(inTen + i) * *(inVec + iVec*nC + i);
                }
                break;
            case 2:
                /* two-fold-loop:
                 *  i - columns
                 *  j - rows
                 */
                for(i = 0; i < nC; i++) {
                for(j = 0; j < nC; j++) {
                    *(outMeas + iVec) = *(outMeas + iVec) + 
                            *(inTen + j + i*nC) * 
                               *(inVec + iVec*nC + i) *
                               *(inVec + iVec*nC + j);
                }}
                break;
            case 3:
                /* three-fold-loop:
                 *  i - pages
                 *  j - columns
                 *  k - rows
                 */
                for(i = 0; i < nC; i++) {
                for(j = 0; j <= i; j++) {
                for(k = 0; k <= j; k++) {
                    *(outMeas + iVec) = *(outMeas + iVec) + 
                            *(inTen + k + j*nC + i*nC*nC) * 
                               *(inVec + iVec*nC + i) *
                               *(inVec + iVec*nC + j) * 
                               *(inVec + iVec*nC + k);
                }}}
                break;
            default:
                return;
        }
    }
    
    return;
}