idamnjanovic@51: /**************************************************************************
idamnjanovic@51:  *
idamnjanovic@51:  * File name: myblas.c
idamnjanovic@51:  *
idamnjanovic@51:  * Ron Rubinstein
idamnjanovic@51:  * Computer Science Department
idamnjanovic@51:  * Technion, Haifa 32000 Israel
idamnjanovic@51:  * ronrubin@cs
idamnjanovic@51:  *
idamnjanovic@51:  * Version: 1.1
idamnjanovic@51:  * Last updated: 13.8.2009
idamnjanovic@51:  *
idamnjanovic@51:  *************************************************************************/
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: #include "myblas.h"
idamnjanovic@51: #include <ctype.h>
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* find maximum of absolute values */
idamnjanovic@51: 
idamnjanovic@51: mwIndex maxabs(double c[], mwSize m)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex maxid=0, k;
idamnjanovic@51:   double absval, maxval = SQR(*c);   /* use square which is quicker than absolute value */
idamnjanovic@51: 
idamnjanovic@51:   for (k=1; k<m; ++k) {
idamnjanovic@51:     absval = SQR(c[k]);
idamnjanovic@51:     if (absval > maxval) {
idamnjanovic@51:       maxval = absval;
idamnjanovic@51:       maxid = k;
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   return maxid;
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* compute y := alpha*x + y */
idamnjanovic@51: 
idamnjanovic@51: void vec_sum(double alpha, double x[], double y[], mwSize n)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex i;
idamnjanovic@51: 
idamnjanovic@51:   for (i=0; i<n; ++i) {
idamnjanovic@51:     y[i] += alpha*x[i];
idamnjanovic@51:   }
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* compute y := alpha*A*x */
idamnjanovic@51: 
idamnjanovic@51: void mat_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex i, j, i_n;
idamnjanovic@51:   double *Ax;
idamnjanovic@51: 
idamnjanovic@51:   Ax = mxCalloc(n,sizeof(double));
idamnjanovic@51: 
idamnjanovic@51:   for (i=0; i<m; ++i) {
idamnjanovic@51:     i_n = i*n;
idamnjanovic@51:     for (j=0; j<n; ++j) {
idamnjanovic@51:       Ax[j] += A[i_n+j] * x[i];
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51: 
idamnjanovic@51:   for (j=0; j<n; ++j) {
idamnjanovic@51:     y[j] = alpha*Ax[j];
idamnjanovic@51:   }
idamnjanovic@51: 
idamnjanovic@51:   mxFree(Ax);
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* compute y := alpha*A'*x */
idamnjanovic@51: 
idamnjanovic@51: void matT_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex i, j, n_i;
idamnjanovic@51:   double sum0, sum1, sum2, sum3;
idamnjanovic@51: 
idamnjanovic@51:   for (j=0; j<m; ++j) {
idamnjanovic@51:     y[j] = 0;
idamnjanovic@51:   }
idamnjanovic@51: 
idamnjanovic@51:   /* use loop unrolling to accelerate computation */
idamnjanovic@51: 
idamnjanovic@51:   for (i=0; i<m; ++i) {
idamnjanovic@51:     n_i = n*i;
idamnjanovic@51:     sum0 = sum1 = sum2 = sum3 = 0;
idamnjanovic@51:     for (j=0; j+4<n; j+=4) {
idamnjanovic@51:       sum0 += A[n_i+j]*x[j];
idamnjanovic@51:       sum1 += A[n_i+j+1]*x[j+1];
idamnjanovic@51:       sum2 += A[n_i+j+2]*x[j+2];
idamnjanovic@51:       sum3 += A[n_i+j+3]*x[j+3];
idamnjanovic@51:     }
idamnjanovic@51:     y[i] += alpha * ((sum0 + sum1) + (sum2 + sum3));
idamnjanovic@51:     while (j<n) {
idamnjanovic@51:       y[i] += alpha*A[n_i+j]*x[j];
idamnjanovic@51:       j++;
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* compute y := alpha*A*x */
idamnjanovic@51: 
idamnjanovic@51: void mat_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m)
idamnjanovic@51: {
idamnjanovic@51:   
idamnjanovic@51:   mwIndex i, j, j1, j2;
idamnjanovic@51: 
idamnjanovic@51:   for (i=0; i<n; ++i) {
idamnjanovic@51:     y[i] = 0;
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   j2 = jc[0];
idamnjanovic@51:   for (i=0; i<m; ++i) {
idamnjanovic@51:     j1 = j2; j2 = jc[i+1];
idamnjanovic@51:     for (j=j1; j<j2; ++j) {
idamnjanovic@51:       y[ir[j]] += alpha * pr[j] * x[i];
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* compute y := alpha*A'*x */
idamnjanovic@51: 
idamnjanovic@51: void matT_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m)
idamnjanovic@51: {
idamnjanovic@51:   
idamnjanovic@51:   mwIndex i, j, j1, j2;
idamnjanovic@51:   
idamnjanovic@51:   for (i=0; i<m; ++i) {
idamnjanovic@51:     y[i] = 0;
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   j2 = jc[0];
idamnjanovic@51:   for (i=0; i<m; ++i) {
idamnjanovic@51:     j1 = j2; j2 = jc[i+1];
idamnjanovic@51:     for (j=j1; j<j2; ++j) {
idamnjanovic@51:       y[i] += alpha * pr[j] * x[ir[j]];
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* compute y := alpha*A*x */
idamnjanovic@51: 
idamnjanovic@51: void mat_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m)
idamnjanovic@51: {
idamnjanovic@51:   
idamnjanovic@51:   mwIndex i, j, j_n, k, kend;
idamnjanovic@51:   
idamnjanovic@51:   for (i=0; i<n; ++i) {
idamnjanovic@51:     y[i] = 0;
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   kend = jc[1];
idamnjanovic@51:   if (kend==0) {   /* x is empty */
idamnjanovic@51:     return;
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   for (k=0; k<kend; ++k) {
idamnjanovic@51:     j = ir[k];
idamnjanovic@51:     j_n = j*n;
idamnjanovic@51:     for (i=0; i<n; ++i) {
idamnjanovic@51:       y[i] += alpha * A[i+j_n] * pr[k];
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51: 
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* compute y := alpha*A'*x */
idamnjanovic@51: 
idamnjanovic@51: void matT_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m)
idamnjanovic@51: {
idamnjanovic@51:   
idamnjanovic@51:   mwIndex i, j, j_n, k, kend;
idamnjanovic@51:   
idamnjanovic@51:   for (i=0; i<m; ++i) {
idamnjanovic@51:     y[i] = 0;
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   kend = jc[1];
idamnjanovic@51:   if (kend==0) {   /* x is empty */
idamnjanovic@51:     return;
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   for (j=0; j<m; ++j) {
idamnjanovic@51:     j_n = j*n;
idamnjanovic@51:     for (k=0; k<kend; ++k) {
idamnjanovic@51:       i = ir[k];
idamnjanovic@51:       y[j] += alpha * A[i+j_n] * pr[k];
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* compute y := alpha*A*x */
idamnjanovic@51: 
idamnjanovic@51: void mat_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m)
idamnjanovic@51: {
idamnjanovic@51:   
idamnjanovic@51:   mwIndex i, j, k, kend, j1, j2;
idamnjanovic@51: 
idamnjanovic@51:   for (i=0; i<n; ++i) {
idamnjanovic@51:     y[i] = 0;
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   kend = jcx[1]; 
idamnjanovic@51:   if (kend==0) {   /* x is empty */
idamnjanovic@51:     return;
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   for (k=0; k<kend; ++k) {
idamnjanovic@51:     i = irx[k];
idamnjanovic@51:     j1 = jc[i]; j2 = jc[i+1];
idamnjanovic@51:     for (j=j1; j<j2; ++j) {
idamnjanovic@51:       y[ir[j]] += alpha * pr[j] * prx[k];
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* compute y := alpha*A'*x */
idamnjanovic@51: 
idamnjanovic@51: void matT_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m)
idamnjanovic@51: {
idamnjanovic@51:   
idamnjanovic@51:   mwIndex i, j, k, jend, kend, jadd, kadd, delta;
idamnjanovic@51:   
idamnjanovic@51:   for (i=0; i<m; ++i) {
idamnjanovic@51:     y[i] = 0;
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   kend = jcx[1];
idamnjanovic@51:   if (kend==0) {   /* x is empty */
idamnjanovic@51:     return;
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   for (i=0; i<m; ++i) {
idamnjanovic@51:     j = jc[i]; 
idamnjanovic@51:     jend = jc[i+1];
idamnjanovic@51:     k = 0;
idamnjanovic@51:     while (j<jend && k<kend) {
idamnjanovic@51:       
idamnjanovic@51:       delta = ir[j] - irx[k];
idamnjanovic@51:       
idamnjanovic@51:       if (delta) { /* if indices differ - increment the smaller one */
idamnjanovic@51:         jadd = delta<0;
idamnjanovic@51:         kadd = 1-jadd;
idamnjanovic@51:         j += jadd;
idamnjanovic@51:         k += kadd;
idamnjanovic@51:       }
idamnjanovic@51:       
idamnjanovic@51:       else {    /* indices are equal - add to result and increment both */
idamnjanovic@51:         y[i] += alpha * pr[j] * prx[k];
idamnjanovic@51:         j++; k++;
idamnjanovic@51:       }
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* matrix-matrix multiplication */
idamnjanovic@51: 
idamnjanovic@51: void mat_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex i1, i2, i3, iX, iA, i2_n;
idamnjanovic@51:   double b;
idamnjanovic@51:   
idamnjanovic@51:   for (i1=0; i1<n*k; i1++) {
idamnjanovic@51:     X[i1] = 0;
idamnjanovic@51:   }
idamnjanovic@51: 
idamnjanovic@51:   for (i2=0; i2<m; ++i2) {
idamnjanovic@51:     i2_n = i2*n;
idamnjanovic@51:     iX = 0;
idamnjanovic@51:     for (i3=0; i3<k; ++i3) {
idamnjanovic@51:       iA = i2_n;
idamnjanovic@51:       b = B[i2+i3*m];
idamnjanovic@51:       for (i1=0; i1<n; ++i1) {
idamnjanovic@51:         X[iX++] += A[iA++]*b;
idamnjanovic@51:       }
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   for (i1=0; i1<n*k; i1++) {
idamnjanovic@51:     X[i1] *= alpha;
idamnjanovic@51:   }
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* matrix-transpose-matrix multiplication */
idamnjanovic@51: 
idamnjanovic@51: void matT_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex i1, i2, i3, iX, iA, i2_n;
idamnjanovic@51:   double *x, sum0, sum1, sum2, sum3;
idamnjanovic@51: 
idamnjanovic@51:   for (i2=0; i2<m; ++i2) {
idamnjanovic@51:     for (i3=0; i3<k; ++i3) {
idamnjanovic@51:       sum0 = sum1 = sum2 = sum3 = 0;
idamnjanovic@51:       for (i1=0; i1+4<n; i1+=4) {
idamnjanovic@51:         sum0 += A[i1+0+i2*n]*B[i1+0+i3*n];
idamnjanovic@51:         sum1 += A[i1+1+i2*n]*B[i1+1+i3*n];
idamnjanovic@51:         sum2 += A[i1+2+i2*n]*B[i1+2+i3*n];
idamnjanovic@51:         sum3 += A[i1+3+i2*n]*B[i1+3+i3*n];
idamnjanovic@51:       }
idamnjanovic@51:       X[i2+i3*m] = (sum0+sum1) + (sum2+sum3);
idamnjanovic@51:       while(i1<n) {
idamnjanovic@51:         X[i2+i3*m] += A[i1+i2*n]*B[i1+i3*n];
idamnjanovic@51:         i1++;
idamnjanovic@51:       }
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   for (i1=0; i1<m*k; i1++) {
idamnjanovic@51:     X[i1] *= alpha;
idamnjanovic@51:   }
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* tensor-matrix product */
idamnjanovic@51: 
idamnjanovic@51: void tens_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex i1, i2, i3, i4, i2_n, nml;
idamnjanovic@51:   double b;
idamnjanovic@51:   
idamnjanovic@51:   nml = n*m*l;
idamnjanovic@51:   for (i1=0; i1<nml; ++i1) {
idamnjanovic@51:     X[i1] = 0;
idamnjanovic@51:   }
idamnjanovic@51: 
idamnjanovic@51:   for (i2=0; i2<m; ++i2) {
idamnjanovic@51:     i2_n = i2*n;
idamnjanovic@51:     for (i3=0; i3<k; ++i3) {
idamnjanovic@51:       for (i4=0; i4<l; ++i4) {
idamnjanovic@51:         b = B[i4+i3*l];
idamnjanovic@51:         for (i1=0; i1<n; ++i1) {
idamnjanovic@51:           X[i1 + i2_n + i4*n*m] += A[i1 + i2_n + i3*n*m] * b;
idamnjanovic@51:         }
idamnjanovic@51:       }
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   for (i1=0; i1<nml; ++i1) {
idamnjanovic@51:     X[i1] *= alpha;
idamnjanovic@51:   }
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* tensor-matrix-transpose product */
idamnjanovic@51: 
idamnjanovic@51: void tens_matT(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex i1, i2, i3, i4, i2_n, nml;
idamnjanovic@51:   double b;
idamnjanovic@51:   
idamnjanovic@51:   nml = n*m*l;
idamnjanovic@51:   for (i1=0; i1<nml; ++i1) {
idamnjanovic@51:     X[i1] = 0;
idamnjanovic@51:   }
idamnjanovic@51: 
idamnjanovic@51:   for (i2=0; i2<m; ++i2) {
idamnjanovic@51:     i2_n = i2*n;
idamnjanovic@51:     for (i4=0; i4<l; ++i4) {
idamnjanovic@51:       for (i3=0; i3<k; ++i3) {
idamnjanovic@51:         b = B[i3+i4*k];
idamnjanovic@51:         for (i1=0; i1<n; ++i1) {
idamnjanovic@51:           X[i1 + i2_n + i4*n*m] += A[i1 + i2_n + i3*n*m] * b;
idamnjanovic@51:         }
idamnjanovic@51:       }
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   for (i1=0; i1<nml; ++i1) {
idamnjanovic@51:     X[i1] *= alpha;
idamnjanovic@51:   }
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* dot product */
idamnjanovic@51: 
idamnjanovic@51: double dotprod(double a[], double b[], mwSize n)
idamnjanovic@51: {
idamnjanovic@51:   double sum = 0;
idamnjanovic@51:   mwIndex i;
idamnjanovic@51:   for (i=0; i<n; ++i)
idamnjanovic@51:     sum += a[i]*b[i];
idamnjanovic@51:   return sum;
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* find maximum of vector */
idamnjanovic@51: 
idamnjanovic@51: mwIndex maxpos(double c[], mwSize m)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex maxid=0, k;
idamnjanovic@51:   double val, maxval = *c;
idamnjanovic@51: 
idamnjanovic@51:   for (k=1; k<m; ++k) {
idamnjanovic@51:     val = c[k];
idamnjanovic@51:     if (val > maxval) {
idamnjanovic@51:       maxval = val;
idamnjanovic@51:       maxid = k;
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   return maxid;
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* solve L*x = b */
idamnjanovic@51: 
idamnjanovic@51: void backsubst_L(double L[], double b[], double x[], mwSize n, mwSize k)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex i, j;
idamnjanovic@51:   double rhs;
idamnjanovic@51: 
idamnjanovic@51:   for (i=0; i<k; ++i) {
idamnjanovic@51:     rhs = b[i];
idamnjanovic@51:     for (j=0; j<i; ++j) {
idamnjanovic@51:       rhs -= L[j*n+i]*x[j];
idamnjanovic@51:     }
idamnjanovic@51:     x[i] = rhs/L[i*n+i];
idamnjanovic@51:   }
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* solve L'*x = b */
idamnjanovic@51: 
idamnjanovic@51: void backsubst_Lt(double L[], double b[], double x[], mwSize n, mwSize k)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex i, j;
idamnjanovic@51:   double rhs;
idamnjanovic@51: 
idamnjanovic@51:   for (i=k; i>=1; --i) {
idamnjanovic@51:     rhs = b[i-1];
idamnjanovic@51:     for (j=i; j<k; ++j) {
idamnjanovic@51:       rhs -= L[(i-1)*n+j]*x[j];
idamnjanovic@51:     }
idamnjanovic@51:     x[i-1] = rhs/L[(i-1)*n+i-1];
idamnjanovic@51:   }
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* solve U*x = b */
idamnjanovic@51: 
idamnjanovic@51: void backsubst_U(double U[], double b[], double x[], mwSize n, mwSize k)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex i, j;
idamnjanovic@51:   double rhs;
idamnjanovic@51: 
idamnjanovic@51:   for (i=k; i>=1; --i) {
idamnjanovic@51:     rhs = b[i-1];
idamnjanovic@51:     for (j=i; j<k; ++j) {
idamnjanovic@51:       rhs -= U[j*n+i-1]*x[j];
idamnjanovic@51:     }
idamnjanovic@51:     x[i-1] = rhs/U[(i-1)*n+i-1];
idamnjanovic@51:   }
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* solve U'*x = b */
idamnjanovic@51: 
idamnjanovic@51: void backsubst_Ut(double U[], double b[], double x[], mwSize n, mwSize k)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex i, j;
idamnjanovic@51:   double rhs;
idamnjanovic@51: 
idamnjanovic@51:   for (i=0; i<k; ++i) {
idamnjanovic@51:     rhs = b[i];
idamnjanovic@51:     for (j=0; j<i; ++j) {
idamnjanovic@51:       rhs -= U[i*n+j]*x[j];
idamnjanovic@51:     }
idamnjanovic@51:     x[i] = rhs/U[i*n+i];
idamnjanovic@51:   }
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* back substitution solver */
idamnjanovic@51: 
idamnjanovic@51: void backsubst(char ul, double A[], double b[], double x[], mwSize n, mwSize k)
idamnjanovic@51: {
idamnjanovic@51:   if (tolower(ul) == 'u') {
idamnjanovic@51:     backsubst_U(A, b, x, n, k);
idamnjanovic@51:   }
idamnjanovic@51:   else if (tolower(ul) == 'l') {
idamnjanovic@51:     backsubst_L(A, b, x, n, k);
idamnjanovic@51:   }
idamnjanovic@51:   else {
idamnjanovic@51:     mexErrMsgTxt("Invalid triangular matrix type: must be ''U'' or ''L''");
idamnjanovic@51:   }
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* solve equation set using cholesky decomposition */
idamnjanovic@51: 
idamnjanovic@51: void cholsolve(char ul, double A[], double b[], double x[], mwSize n, mwSize k)
idamnjanovic@51: {
idamnjanovic@51:   double *tmp;
idamnjanovic@51: 
idamnjanovic@51:   tmp = mxMalloc(k*sizeof(double));
idamnjanovic@51: 
idamnjanovic@51:   if (tolower(ul) == 'l') {
idamnjanovic@51:     backsubst_L(A, b, tmp, n, k);
idamnjanovic@51:     backsubst_Lt(A, tmp, x, n, k);
idamnjanovic@51:   }
idamnjanovic@51:   else if (tolower(ul) == 'u') {
idamnjanovic@51:     backsubst_Ut(A, b, tmp, n, k);
idamnjanovic@51:     backsubst_U(A, tmp, x, n, k);
idamnjanovic@51:   }
idamnjanovic@51:   else {
idamnjanovic@51:     mexErrMsgTxt("Invalid triangular matrix type: must be either ''U'' or ''L''");
idamnjanovic@51:   }
idamnjanovic@51: 
idamnjanovic@51:   mxFree(tmp);
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* perform a permutation assignment y := x(ind(1:k)) */
idamnjanovic@51: 
idamnjanovic@51: void vec_assign(double y[], double x[], mwIndex ind[], mwSize k)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex i;
idamnjanovic@51: 
idamnjanovic@51:   for (i=0; i<k; ++i)
idamnjanovic@51:     y[i] = x[ind[i]];
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* matrix transpose */
idamnjanovic@51: 
idamnjanovic@51: void transpose(double X[], double Y[], mwSize n, mwSize m)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex i, j, i_m, j_n;
idamnjanovic@51:   
idamnjanovic@51:   if (n<m) {
idamnjanovic@51:     for (j=0; j<m; ++j) {
idamnjanovic@51:       j_n = j*n;
idamnjanovic@51:       for (i=0; i<n; ++i) {
idamnjanovic@51:         Y[j+i*m] = X[i+j_n];
idamnjanovic@51:       }
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   else {
idamnjanovic@51:     for (i=0; i<n; ++i) {
idamnjanovic@51:       i_m = i*m;
idamnjanovic@51:       for (j=0; j<m; ++j) {
idamnjanovic@51:         Y[j+i_m] = X[i+j*n];
idamnjanovic@51:       }
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* print contents of matrix */
idamnjanovic@51: 
idamnjanovic@51: void printmat(double A[], int n, int m, char* matname)
idamnjanovic@51: {
idamnjanovic@51:   int i, j;
idamnjanovic@51:   mexPrintf("\n%s = \n\n", matname);
idamnjanovic@51: 
idamnjanovic@51:   if (n*m==0) {
idamnjanovic@51:     mexPrintf("   Empty matrix: %d-by-%d\n\n", n, m);
idamnjanovic@51:     return;
idamnjanovic@51:   }
idamnjanovic@51: 
idamnjanovic@51:   for (i=0; i<n; ++i) {
idamnjanovic@51:     for (j=0; j<m; ++j)
idamnjanovic@51:       mexPrintf("   %lf", A[j*n+i]);
idamnjanovic@51:     mexPrintf("\n");
idamnjanovic@51:   }
idamnjanovic@51:   mexPrintf("\n");
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* print contents of sparse matrix */
idamnjanovic@51: 
idamnjanovic@51: void printspmat(mxArray *a, char* matname)
idamnjanovic@51: {
idamnjanovic@51:   mwIndex *aJc = mxGetJc(a);
idamnjanovic@51:   mwIndex *aIr = mxGetIr(a);
idamnjanovic@51:   double *aPr = mxGetPr(a);
idamnjanovic@51: 
idamnjanovic@51:   int i;
idamnjanovic@51: 
idamnjanovic@51:   mexPrintf("\n%s = \n\n", matname);
idamnjanovic@51: 
idamnjanovic@51:   for (i=0; i<aJc[1]; ++i)
idamnjanovic@51:     printf("   (%d,1) = %lf\n", aIr[i]+1,aPr[i]);
idamnjanovic@51: 
idamnjanovic@51:   mexPrintf("\n");
idamnjanovic@51: }
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: /* matrix multiplication using Winograd's algorithm */
idamnjanovic@51: 
idamnjanovic@51: /*
idamnjanovic@51: void mat_mat2(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k)
idamnjanovic@51: {
idamnjanovic@51:   
idamnjanovic@51:   mwIndex i1, i2, i3, iX, iA, i2_n;
idamnjanovic@51:   double b, *AA, *BB;
idamnjanovic@51:   
idamnjanovic@51:   AA = mxCalloc(n,sizeof(double));
idamnjanovic@51:   BB = mxCalloc(k,sizeof(double));
idamnjanovic@51:   
idamnjanovic@51:   for (i1=0; i1<n*k; i1++) {
idamnjanovic@51:     X[i1] = 0;
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   for (i1=0; i1<n; ++i1) {
idamnjanovic@51:     for (i2=0; i2<m/2; ++i2) {
idamnjanovic@51:       AA[i1] += A[i1+2*i2*n]*A[i1+(2*i2+1)*n];
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51: 
idamnjanovic@51:   for (i2=0; i2<k; ++i2) {
idamnjanovic@51:     for (i1=0; i1<m/2; ++i1) {
idamnjanovic@51:       BB[i2] += B[2*i1+i2*m]*B[2*i1+1+i2*m];
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51: 
idamnjanovic@51:   for (i2=0; i2<k; ++i2) {
idamnjanovic@51:     for (i3=0; i3<m/2; ++i3) {
idamnjanovic@51:       for (i1=0; i1<n; ++i1) {
idamnjanovic@51:         X[i1+i2*n] += (A[i1+(2*i3)*n]+B[2*i3+1+i2*m])*(A[i1+(2*i3+1)*n]+B[2*i3+i2*m]);
idamnjanovic@51:       }
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   if (m%2) {
idamnjanovic@51:     for (i2=0; i2<k; ++i2) {
idamnjanovic@51:       for (i1=0; i1<n; ++i1) {
idamnjanovic@51:         X[i1+i2*n] += A[i1+(m-1)*n]*B[m-1+i2*m];
idamnjanovic@51:       }
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   for (i2=0; i2<k; ++i2) {
idamnjanovic@51:     for (i1=0; i1<n; ++i1) {
idamnjanovic@51:       X[i1+i2*n] -= (AA[i1] + BB[i2]);
idamnjanovic@51:       X[i1+i2*n] *= alpha;
idamnjanovic@51:     }
idamnjanovic@51:   }
idamnjanovic@51:   
idamnjanovic@51:   mxFree(AA);
idamnjanovic@51:   mxFree(BB);
idamnjanovic@51: }
idamnjanovic@51: */
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: 
idamnjanovic@51: