text:
// lu2tlub.c
//
//  blocked LU decomposition library for column-major version
//
// Time-stamp: <11/05/13 11:35:56 makino>

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifndef NOBLAS
#ifdef MKL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#endif

#include <lu2lib.h>

typedef double v2df __attribute__((vector_size(16)));
typedef union {v2df v; double s[2];}v2u;

void timer_init();

#define MAXTHREADS 4

static int findpivot0(int n, double a[][n], int current)
{
    double amax = fabs(a[0][current]);
    int i;
    int p=current;
    BEGIN_TSC;
    for(i=current+1;i<n;i++){
	if (fabs(a[0][i]) > amax){
	    amax = fabs(a[0][i]);
	    p = i;
	}
    }
    END_TSC(t,7);
    return p;
}
static int findpivot(int n, double a[][n], int current)
{
    int p;
    int p2;
    BEGIN_TSC;
    p = cblas_idamax(n-current, a[0]+current, 1)+current;
    //    p2 = findpivot0( n, a, current);
    //    printf("n, current, p, p2 = %d %d %d %d\n",n, current,p, p2);
    END_TSC(t,7);
    return p;
}
static int findpivot_sequentical(int n, double a[][n], int current)
{
    double amax = fabs(a[0][current]);
    int i;
    int p=current;
    for(i=current+1;i<n;i++){
	if (fabs(a[0][i]) > amax){
	    amax = fabs(a[0][i]);
	    p = i;
	}
    }
    return p;
}
static int findpivot_omp(int n, double a[][n], int current)
// factor 2 slower than sequential code even for n=8k....
// on Core i7 920
{
    double amax[MAXTHREADS]={-1.0,-1.0,-1.0,-1.0};
    int p[MAXTHREADS];
    double am;
    int pm;
    int di= (n-current-1)/4;
    int k;
    BEGIN_TSC;
    if (di > 1024){
	//#pragma omp parallel for private(k)
	for (k=0;k<MAXTHREADS; k++){
	    int i;
	    int istart = current+1+k*di;
	    int iend = current+1+(k+1)*di;
	    if (iend > n) iend = n;
	    for(i=istart;i<iend;i++){
		if (fabs(a[0][i]) > amax[k]){
		    amax[k] = fabs(a[0][i]);
		    p[k] = i;
		}
	    }
	}
	pm =p[0];
	am = amax[0];
	for (k=1;k<MAXTHREADS; k++){
	    if(amax[k]>am){
		am = amax[k];
		pm=p[k];
	    }
	}
    }else{
	pm=findpivot_sequentical( n,  a, current);
    }
    END_TSC(t,7);
    return pm;
}

static void swaprows(int n, double a[][n], int row1, int row2,
		     int cstart, int cend)
{
    int j;
    if (row1 != row2){
	for(j=cstart;j<cend;j++){
	    double tmp = a[j][row1];
	    a[j][row1] = a[j][row2];
	    a[j][row2]=tmp;
	}
    }
}

static void scalerow( int n, double a[n+1][n], double scale,
	       int row, int cstart, int cend)
{
    int j;
    for(j=cstart;j<cend;j++) a[j][row]*= scale;
}

static void vsmulandsub(int n, double a[n+1][n], int cr, int cc,
		 int c0,  int r0,int r1)
{
    int j,k;
    double * ar = a[cr];
    k=c0;
    double s = a[c0][cc];
    double *al = a[c0];
    while (r0 & 7){
	al[r0] -= ar[r0]*s;
	r0++;
    }
    while (r1 & 7){
	al[r1-1] -= ar[r1-1]*s;
	r1--;
    }
    v2df * arv = (v2df*) (ar+r0);
    v2df * alv = (v2df*) (al+r0);
    v2df ss = (v2df){s,s};
    //    for(j=r0;j<r1;j++)
    //	al[j] -= ar[j]*s;
    for(j=0;j<(r1-r0)/2;j+=4){
	alv[j] -= arv[j]*ss;
	alv[j+1] -= arv[j+1]*ss;
	alv[j+2] -= arv[j+2]*ss;
	alv[j+3] -= arv[j+3]*ss;
	__builtin_prefetch(alv+j+32,1,3);
	__builtin_prefetch(arv+j+32,0,0);

    }
}

static void vvmulandsub(int n, double a[n+1][n], int cr, int cc,
		 int c0, int c1, int r0,int r1)
{
    int j,k;
    double * ar = a[cr];
#ifdef TIMETEST
    BEGIN_TSC;
#endif
    if (c1-c0 == 1){
	vsmulandsub(n,a, cr, cc,c0, r0,r1);
    }else{
	for (k=c0;k<c1;k++){
	    double s = a[k][cc];
	    double *al = a[k];
	    for(j=r0;j<r1;j++)
		al[j] -= ar[j]*s;
	}
    }
#ifdef TIMETEST
    END_TSC(t,6);
#endif
}
static void mmmulandsub_old(int n, double a[n+1][n], int m0, int m1,
		 int c0, int c1, int r0,int r1)
{
    int j,k,l;
    printf("Enter mmul\n");
#ifndef NOBLAS
    cblas_dgemm( CblasColMajor, CblasNoTrans, CblasNoTrans,
		 r1-r0, c1-c0, m1-m0, -1.0, &(a[m0][r0]), n,
		 &(a[c0][m0]), n, 1, &(a[c0][r0]), n );
    // example:
    // r0, m0 = i+m,i
    // m0, c0 = i, i+m
    // r0, c0 = i+m, i+m
    //r1-r0 = n-i-m
    // c1-c0 = iend-i-m
    // m1-m0 = m
#else
    for(j=r0;j<r1;j++)
	for (k=c0;k<c1;k++)
	    for (l=m0; l<m1; l++)
		a[k][j] -= a[l][j]*a[k][l];
#endif
}

static void matmul_for_nk4_0(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int j,k,l;
    for(l=0;l<4;l+=2){
	for(k=0;k<4;k+=2){
	    for(j=0;j<n;j++){
		c[k][j] -= a[l][j]*b[k][l]+ a[l+1][j]*b[k][l+1];
		c[k+1][j] -= a[l][j]*b[k+1][l]+ a[l+1][j]*b[k+1][l+1];
	    }
	}
    }
}
static void matmul_for_nk8_0(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int i,j,k,l;
    const int m=8;
    double btmp[m][m];

    for(i=0;i<m;i++)for(j=0;j<m;j++)btmp[i][j]=b[i][j];
#pragma omp parallel for private(l,k,j)
    for(l=0;l<8;l+=2){
	for(k=0;k<8;k+=2){
	    for(j=0;j<n;j++){
		c[k][j] -= a[l][j]*btmp[k][l]+ a[l+1][j]*btmp[k][l+1];
		c[k+1][j] -= a[l][j]*btmp[k+1][l]+ a[l+1][j]*btmp[k+1][l+1];
	    }
	}
    }

}
static void matmul_for_nk8_2(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int i,j,k,l;
    const int m=8;
    double btmp[m][m];

    for(i=0;i<m;i++)for(j=0;j<m;j++)btmp[i][j]=b[i][j];
#pragma omp parallel for private(l,k,j)
    for(k=0;k<8;k+=2){
	for(l=0;l<8;l+=4){
	    for(j=0;j<n;j++){
		c[k][j] -= a[l][j]*btmp[k][l]
		    + a[l+1][j]*btmp[k][l+1]
		    + a[l+2][j]*btmp[k][l+2]
		    + a[l+3][j]*btmp[k][l+3];
		c[k+1][j] -= a[l][j]*btmp[k+1][l]
		    + a[l+1][j]*btmp[k+1][l+1]
		    + a[l+2][j]*btmp[k+1][l+2]
		    + a[l+3][j]*btmp[k+1][l+3];
		//		c[k][j] -= a[l][j]*btmp[k][l]+ a[l+1][j]*btmp[k][l+1];
		//		c[k+1][j] -= a[l][j]*btmp[k+1][l]+ a[l+1][j]*btmp[k+1][l+1];
	    }
	}
    }

}
static void matmul_for_nk8_3(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int i,j,k,l;
    const int m=8;
    const int mm = 32;
    double btmp[m][m];
    double atmp[m][mm];

    for(i=0;i<m;i++)for(j=0;j<m;j++)btmp[i][j]=b[i][j];
    for(j=0;j<n;j+=mm){
	int jj;
	for(i=0;i<m;i++){
	    for(k=0;k<mm;k++){
		atmp[i][k]=a[i][j+k];
	    }
	}
	for(i=0;i<m;i++){
	    int jjend = mm;
	    if (jjend+j > n) jjend = n-j;
	    for(jj=0;jj<jjend;jj+=2){
		c[i][j+jj] -= atmp[0][jj]*btmp[i][0]
		    +atmp[1][jj]*btmp[i][1]
		    +atmp[2][jj]*btmp[i][2]
		    +atmp[3][jj]*btmp[i][3]
		    +atmp[4][jj]*btmp[i][4]
		    +atmp[5][jj]*btmp[i][5]
		    +atmp[6][jj]*btmp[i][6]
		    +atmp[7][jj]*btmp[i][7];
		c[i][j+jj+1] -= atmp[0][jj+1]*btmp[i][0]
		    +atmp[1][jj+1]*btmp[i][1]
		    +atmp[2][jj+1]*btmp[i][2]
		    +atmp[3][jj+1]*btmp[i][3]
		    +atmp[4][jj+1]*btmp[i][4]
		    +atmp[5][jj+1]*btmp[i][5]
		    +atmp[6][jj+1]*btmp[i][6]
		    +atmp[7][jj+1]*btmp[i][7];
	    }
	}
    }

}
static void matmul_for_nk8_4(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int i,j,k,l;
    const int m=8;
    const int mm = 16;
    double btmp[m][m];
    v2df b2tmp[m][m];
    double atmp[m][mm];

    for(i=0;i<m;i++)for(j=0;j<m;j++)b2tmp[i][j]=(v2df){b[i][j],b[i][j]};
    //#pragma omp parallel for private(j,i,atmp)
    //use of OMP here does not speed things up
    for(j=0;j<n;j+=mm){
	int jj;
	for(i=0;i<m;i++){
	    v2df * dest = (v2df*)(atmp[i]);
	    v2df * src = (v2df*)(a[i]+j);
	    dest[0]=src[0];
	    dest[ 1]=src[ 1];
	    dest[ 2]=src[ 2];
	    dest[ 3]=src[ 3];
	    dest[ 4]=src[ 4];
	    dest[ 5]=src[ 5];
	    dest[ 6]=src[ 6];
	    dest[ 7]=src[ 7];
	}
	for(i=0;i<m;i++){
	    int jjend = mm;
	    __builtin_prefetch(c[i]+j+64,1,0);
	    __builtin_prefetch(c[i]+j+80,1,0);
	    if (jjend+j > n) jjend = n-j;
	    for(jj=0;jj<jjend;jj+=4){
		v2df* cp = (v2df*)(&c[i][j+jj]);
		v2df* ap = (v2df*)(atmp[0]+jj);
		v2df* cpp = (v2df*)(&c[i][j+jj+2]);
		v2df* app = (v2df*)(atmp[0]+jj+2);
		*cp -= (*ap)*b2tmp[i][0]
		    +(*(ap+16))*b2tmp[i][1]
		    +(*(ap+32))*b2tmp[i][2]
		    +(*(ap+48))*b2tmp[i][3]
		    +(*(ap+64))*b2tmp[i][4]
		    +(*(ap+80))*b2tmp[i][5]
		    +(*(ap+96))*b2tmp[i][6]
		    +(*(ap+112))*b2tmp[i][7];
		*cpp -= (*app)*b2tmp[i][0]
		    +(*(app+16))*b2tmp[i][1]
		    +(*(app+32))*b2tmp[i][2]
		    +(*(app+48))*b2tmp[i][3]
		    +(*(app+64))*b2tmp[i][4]
		    +(*(app+80))*b2tmp[i][5]
		    +(*(app+96))*b2tmp[i][6]
		    +(*(app+112))*b2tmp[i][7];
	    }
	}
    }

}
static void matmul_for_nk8_5(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int i,j,k,l;
    const int m=8;
    const int mm = 32;
    double btmp[m][m];
    v2df b2tmp[m][m];

    for(i=0;i<m;i++)for(j=0;j<m;j++)b2tmp[i][j]=(v2df){b[i][j],b[i][j]};
    //#pragma omp parallel for private(j,i)
    //use of OMP here does not speed things up
    for(j=0;j<n;j+=mm){
	double atmp[m][mm];
	int jj;
	for(i=0;i<m;i++){
	    v2df * dest = (v2df*)(atmp[i]);
	    v2df * src = (v2df*)(a[i]+j);
	    dest[0]=src[0];
	    dest[ 1]=src[ 1];
	    dest[ 2]=src[ 2];
	    dest[ 3]=src[ 3];
	    dest[ 4]=src[ 4];
	    dest[ 5]=src[ 5];
	    dest[ 6]=src[ 6];
	    dest[ 7]=src[ 7];
	    dest[ 8]=src[ 8];
	    dest[ 9]=src[ 9];
	    dest[ 10]=src[ 10];
	    dest[ 11]=src[ 11];
	    dest[ 12]=src[ 12];
	    dest[ 13]=src[ 13];
	    dest[ 14]=src[ 14];
	    dest[ 15]=src[ 15];
	}
	for(i=0;i<m;i++){
	    int jjend = mm;
	    __builtin_prefetch(c[i]+j+64,1,0);
	    __builtin_prefetch(c[i]+j+80,1,0);
	    if (jjend+j > n) jjend = n-j;
	    for(jj=0;jj<jjend;jj+=4){
		v2df* cp = (v2df*)(&c[i][j+jj]);
		v2df* ap = (v2df*)(atmp[0]+jj);
		v2df* cpp = (v2df*)(&c[i][j+jj+2]);
		v2df* app = (v2df*)(atmp[0]+jj+2);
		*cp -= (*ap)*b2tmp[i][0]
		    +(*(ap+16))*b2tmp[i][1]
		    +(*(ap+32))*b2tmp[i][2]
		    +(*(ap+48))*b2tmp[i][3]
		    +(*(ap+64))*b2tmp[i][4]
		    +(*(ap+80))*b2tmp[i][5]
		    +(*(ap+96))*b2tmp[i][6]
		    +(*(ap+112))*b2tmp[i][7];
		*cpp -= (*app)*b2tmp[i][0]
		    +(*(app+16))*b2tmp[i][1]
		    +(*(app+32))*b2tmp[i][2]
		    +(*(app+48))*b2tmp[i][3]
		    +(*(app+64))*b2tmp[i][4]
		    +(*(app+80))*b2tmp[i][5]
		    +(*(app+96))*b2tmp[i][6]
		    +(*(app+112))*b2tmp[i][7];
	    }
	}
    }

}

static void matmul_for_nk8(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int m0,m1,m2,m3;
    int s1, s2, s3;
    int ds = (n/64)*16;
    s1 = ds;
    s2 = ds*2;
    s3 = ds*3;
    m0 = ds;
    m1 = ds;
    m2 = ds;
    m3 = n-s3;
    //    fprintf(stderr,"n, s, m = %d %d %d %d %d %d %d %d\n",
    //    n,s1,s2,s3,m0,m1,m2,m3);

#pragma omp parallel
#pragma omp sections
  {
#pragma omp section
      matmul_for_nk8_5(n1,a,n2,b,n3,c,m0);
#pragma omp section
      matmul_for_nk8_5(n1,((double*)a)+s1 ,n2,b,n3,((double*)c)+s1,m1);
#pragma omp section
      matmul_for_nk8_5(n1,((double*)a)+s2 ,n2,b,n3,((double*)c)+s2,m2);
#pragma omp section
      matmul_for_nk8_5(n1,((double*)a)+s3 ,n2,b,n3,((double*)c)+s3,m3);
  }
}

static void matmul_for_nk4_3(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int i,j,k,l;
    const int m=4;
    const int mm = 16;
    double btmp[m][m];
    v2df b2tmp[m][m];
    double atmp[m][mm];
    v2df atmpt[mm/2][m];

    for(i=0;i<m;i++)for(j=0;j<m;j++)b2tmp[i][j]=(v2df){b[i][j],b[i][j]};
    //#pragma omp parallel for private(j,i,atmp)
    //use of OMP here does not speed things up
    for(j=0;j<n;j+=mm){
	int jj;

	for(i=0;i<m;i++){
	    v2df * dest = (v2df*)(atmp[i]);
	    v2df * src = (v2df*)(a[i]+j);
	    atmpt[ 0][i]=src[0];
	    atmpt[ 1][i]=src[ 1];
	    atmpt[ 2][i]=src[ 2];
	    atmpt[ 3][i]=src[ 3];
	    atmpt[ 4][i]=src[ 4];
	    atmpt[ 5][i]=src[ 5];
	    atmpt[ 6][i]=src[ 6];
	    atmpt[ 7][i]=src[ 7];
	}
	for(i=0;i<m;i++){
	    int jjend = mm;
	    if (jjend+j > n) jjend = n-j;
	    v2df* cp = (v2df*)(&c[i][j]);
	    v2df* ap = (v2df*)(atmpt[0]);
	    v2df* cpp = (v2df*)(&c[i][j+2]);
	    v2df* app = (v2df*)(atmpt[1]);
	    v2df* bp = b2tmp[i];
	    for(jj=0;jj<jjend;jj+=4){
		*cp -= ap[0]*bp[0]
		    +ap[1]*bp[1]
		    +ap[2]*bp[2]
		    +ap[3]*bp[3];
		*cpp -= app[0]*bp[0]
		    +app[1]*bp[1]
		    +app[2]*bp[2]
		    +app[3]*bp[3];
		cp += 2;
		cpp+= 2;
		ap += m*2;
		app+= m*2;
	    }
	}

    }

}

static void matmul_for_nk4_4(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int i,j,k,l;
    const int m=4;
    const int mm = 32;
    double btmp[m][m];
    v2df b2tmp[m][m];
    double atmp[m][mm];

    for(i=0;i<m;i++)for(j=0;j<m;j++)b2tmp[i][j]=(v2df){b[i][j],b[i][j]};
    //#pragma omp parallel for private(j,i,atmp)
    // use of OMP does not speed things up here...
    for(j=0;j<n;j+=mm){
	int jj;
	for(i=0;i<m;i++){
	    v2df * dest = (v2df*)(atmp[i]);
	    v2df * src = (v2df*)(a[i]+j);
#if 0
	    for(k=0;k<mm;k++)atmp[i][k]=a[i][j+k];
#endif
#if 0
	    for(k=0;k<mm/2;k++)dest[k]=src[k];
#endif
	    for(k=0;k<mm/2;k+=2){
		dest[k]=src[k];
		dest[k+1]=src[k+1];
	    }
	}
	for(i=0;i<m;i++){
	    int jjend = mm;
	    v2df * bp = b2tmp[i];
	    if (jjend+j > n) jjend = n-j;
	    for(jj=0;jj<jjend;jj+=2){
		v2df* cp = (v2df*)(&c[i][j+jj]);
		v2df* ap = (v2df*)(atmp[0]+jj);
		v2df* cpp = (v2df*)(&c[i][j+jj+2]);
		v2df* app = (v2df*)(atmp[0]+jj+2);
		*cp -= (*ap)*bp[0]
		    +(*(ap+16))*bp[1]
		    +(*(ap+32))*bp[2]
		    +(*(ap+48))*bp[3];
#if 0
		*cpp -= (*app)*bp[0]
		    +(*(app+16))*bp[1]
		    +(*(app+32))*bp[2]
		    +(*(app+48))*bp[3];
#endif
	    }
	}
    }

}
static void matmul_for_nk4_1(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int i,j,k,l;
    const int m=4;
    const int mm = 32;
    double btmp[m][m];
    v2df b2tmp[m][m];

    for(i=0;i<m;i++)for(j=0;j<m;j++)b2tmp[i][j]=(v2df){b[i][j],b[i][j]};
#pragma omp parallel for private(j,i)
    // use of OMP does not speed things up here...
    for(j=0;j<n;j+=mm){
	double atmp[m][mm];
	int jj;
	for(i=0;i<m;i++){
	    v2df * dest = (v2df*)(atmp[i]);
	    v2df * src = (v2df*)(a[i]+j);
#if 0
	    for(k=0;k<mm;k++)atmp[i][k]=a[i][j+k];
#endif
#if 0
	    for(k=0;k<mm/2;k++)dest[k]=src[k];
#endif
	    for(k=0;k<mm/2;k+=4){
		dest[k]=src[k];
		dest[k+1]=src[k+1];
		dest[k+2]=src[k+2];
		dest[k+3]=src[k+3];
	    }
	}
	for(i=0;i<m;i++){
	    int jjend = mm;
	    v2df * bp = b2tmp[i];
	    if (jjend+j > n) jjend = n-j;
	    for(jj=0;jj<jjend;jj+=4){
		v2df* cp = (v2df*)(&c[i][j+jj]);
		v2df* ap = (v2df*)(atmp[0]+jj);
		v2df* cpp = (v2df*)(&c[i][j+jj+2]);
		v2df* app = (v2df*)(atmp[0]+jj+2);
		*cp -= (*ap)*bp[0]
		    +(*(ap+16))*bp[1]
		    +(*(ap+32))*bp[2]
		    +(*(ap+48))*bp[3];
		*cpp -= (*app)*bp[0]
		    +(*(app+16))*bp[1]
		    +(*(app+32))*bp[2]
		    +(*(app+48))*bp[3];
	    }
	}
    }

}
static void matmul_for_nk2_0(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int j,k,l;
    for(j=0;j<n;j++){
	    c[0][j] -= a[0][j]*b[0][0]+a[1][j]*b[0][1];
	    c[1][j] -= a[0][j]*b[1][0]+a[1][j]*b[1][1];
    }
}
static void matmul_for_nk2_1(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int j,k,l;

    for(j=0;j<n;j+=2){
	    c[0][j] -= a[0][j]*b[0][0]+a[1][j]*b[0][1];
	    c[1][j] -= a[0][j]*b[1][0]+a[1][j]*b[1][1];
	    c[0][j+1] -= a[0][j+1]*b[0][0]+a[1][j+1]*b[0][1];
	    c[1][j+1] -= a[0][j+1]*b[1][0]+a[1][j+1]*b[1][1];
    }
}
static void matmul_for_nk2_2(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int j,k,l;

    register double b00 = b[0][0];
    register double b01 = b[0][1];
    register double b10 = b[1][0];
    register double b11 = b[1][1];
    for(j=0;j<n;j+=2){
	    c[0][j] -= a[0][j]*b00+a[1][j]*b01;
	    c[0][j+1] -= a[0][j+1]*b00+a[1][j+1]*b01;
	    c[1][j] -= a[0][j]*b10+a[1][j]*b11;
	    c[1][j+1] -= a[0][j+1]*b10+a[1][j+1]*b11;
    }

}
static void matmul_for_nk2(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int j,k,l;
    register v2df b00 = (v2df){b[0][0],b[0][0]};
    register v2df b01 = (v2df){b[0][1],b[0][1]};
    register v2df b10 = (v2df){b[1][0],b[1][0]};
    register v2df b11 = (v2df){b[1][1],b[1][1]};
    v2df * a0 = (v2df*) a[0];
    v2df * a1 = (v2df*) a[1];
    v2df * c0 = (v2df*) c[0];
    v2df * c1 = (v2df*) c[1];
    int nh = n>>1;
    if (nh & 1){
	j=nh-1;
	c0[j] -= a0[j]*b00+a1[j]*b01;
	c1[j] -= a0[j]*b10+a1[j]*b11;
	nh = nh-1;
    }
    for(j=0;j<nh;j+=2){
	c0[j] -= a0[j]*b00+a1[j]*b01;
	c0[j+1] -= a0[j+1]*b00+a1[j+1]*b01;
	c1[j] -= a0[j]*b10+a1[j]*b11;
	c1[j+1] -= a0[j+1]*b10+a1[j+1]*b11;
	//	__builtin_prefetch((double*)&a0[j+32],0);
	//	__builtin_prefetch((double*)&a1[j+32],0);
	//	__builtin_prefetch((double*)&c0[j+32],1);
	//	__builtin_prefetch((double*)&c1[j+32],1);
	//	asm("prefetcht2 %0"::"m"(a0[j+32]):"memory");
	//	asm("prefetcht2 %0"::"m"(a1[j+32]):"memory");
    }
}

static void matmul_for_nk4_5(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{

    matmul_for_nk2(n1, a, n2, b, n3,c, n);
    matmul_for_nk2(n1, (double(*)[])  (a[2]), n2,
		   (double(*)[])  &(b[0][2]), n3,c, n);
    matmul_for_nk2(n1, a, n2, (double(*)[])  b[2], n3,
		   (double(*)[])  c[2], n);
    matmul_for_nk2(n1, (double(*)[])  &(a[2]), n2,
		   (double(*)[])  &(b[2][2]), n3,
		   (double(*)[])  c[2], n);

}

static void matmul_for_nk4_6(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    register v2df b00 = (v2df){b[0][0],b[0][0]};
    register v2df b01 = (v2df){b[0][1],b[0][1]};
    register v2df b02 = (v2df){b[0][2],b[0][2]};
    register v2df b03 = (v2df){b[0][3],b[0][3]};
    register v2df b10 = (v2df){b[1][0],b[1][0]};
    register v2df b11 = (v2df){b[1][1],b[1][1]};
    register v2df b12 = (v2df){b[1][2],b[1][2]};
    register v2df b13 = (v2df){b[1][3],b[1][3]};
    register v2df b20 = (v2df){b[2][0],b[2][0]};
    register v2df b21 = (v2df){b[2][1],b[2][1]};
    register v2df b22 = (v2df){b[2][2],b[2][2]};
    register v2df b23 = (v2df){b[2][3],b[2][3]};
    register v2df b30 = (v2df){b[3][0],b[3][0]};
    register v2df b31 = (v2df){b[3][1],b[3][1]};
    register v2df b32 = (v2df){b[3][2],b[3][2]};
    register v2df b33 = (v2df){b[3][3],b[3][3]};

    v2df * a0 = (v2df*) a[0];
    v2df * a1 = (v2df*) a[1];
    v2df * a2 = (v2df*) a[2];
    v2df * a3 = (v2df*) a[3];
    v2df * c0 = (v2df*) c[0];
    v2df * c1 = (v2df*) c[1];
    v2df * c2 = (v2df*) c[2];
    v2df * c3 = (v2df*) c[3];
    int nh = n>>1;
    //#pragma omp parallel
    {
	//#pragma omp section
	{
	    int j;
	    //#pragma omp for private (j)
	    for(j=0;j<nh;j++){
		c0[j] -= a0[j]*b00+a1[j]*b01+a2[j]*b02+a3[j]*b03;
		c1[j] -= a0[j]*b10+a1[j]*b11+a2[j]*b12+a3[j]*b13;
	    }
	}
	//#pragma omp section
	{
	    int j;
	    //#pragma omp for private (j)
	    for(j=0;j<nh;j++){
		c2[j] -= a0[j]*b20+a1[j]*b21+a2[j]*b22+a3[j]*b23;
		c3[j] -= a0[j]*b30+a1[j]*b31+a2[j]*b32+a3[j]*b33;
	    }
	}
    }
}
static void matmul_for_nk4_7(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    v2df * a0 = (v2df*) a[0];
    v2df * a1 = (v2df*) a[1];
    v2df * a2 = (v2df*) a[2];
    v2df * a3 = (v2df*) a[3];
    int j;
    {
	int nh = n>>1;
	v2df * c0 = (v2df*) c[0];
	v2df * c1 = (v2df*) c[1];
	 v2df b00 = (v2df){b[0][0],b[0][0]};
	 v2df b01 = (v2df){b[0][1],b[0][1]};
	 v2df b02 = (v2df){b[0][2],b[0][2]};
	 v2df b03 = (v2df){b[0][3],b[0][3]};
	 v2df b10 = (v2df){b[1][0],b[1][0]};
	 v2df b11 = (v2df){b[1][1],b[1][1]};
	 v2df b12 = (v2df){b[1][2],b[1][2]};
	 v2df b13 = (v2df){b[1][3],b[1][3]};
	if (nh & 1){
	    j=nh-1;
	    c0[j] -= a0[j]*b00+a1[j]*b01+a2[j]*b02+a3[j]*b03;
	    c1[j] -= a0[j]*b10+a1[j]*b11+a2[j]*b12+a3[j]*b13;
	    nh--;
	}
	for(j=0;j<nh;j+=2){
	    c0[j] -= a0[j]*b00+a1[j]*b01+a2[j]*b02+a3[j]*b03;
	    c1[j] -= a0[j]*b10+a1[j]*b11+a2[j]*b12+a3[j]*b13;
	    c0[j+1] -= a0[j+1]*b00+a1[j+1]*b01+a2[j+1]*b02+a3[j+1]*b03;
	    c1[j+1] -= a0[j+1]*b10+a1[j+1]*b11+a2[j+1]*b12+a3[j+1]*b13;
	}
    }
    {
	int nh = n>>1;
	v2df * c2 = (v2df*) c[2];
	v2df * c3 = (v2df*) c[3];
	 v2df b20 = (v2df){b[2][0],b[2][0]};
	 v2df b21 = (v2df){b[2][1],b[2][1]};
	 v2df b22 = (v2df){b[2][2],b[2][2]};
	 v2df b23 = (v2df){b[2][3],b[2][3]};
	register v2df b30 = (v2df){b[3][0],b[3][0]};
	register v2df b31 = (v2df){b[3][1],b[3][1]};
	register v2df b32 = (v2df){b[3][2],b[3][2]};
	register v2df b33 = (v2df){b[3][3],b[3][3]};

	for(j=0;j<nh;j++){
	    c2[j] -= a0[j]*b20+a1[j]*b21+a2[j]*b22+a3[j]*b23;
	    c3[j] -= a0[j]*b30+a1[j]*b31+a2[j]*b32+a3[j]*b33;
	}
    }
}

static void matmul_for_nk4_8(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    v2df * a0 = (v2df*) a[0];
    v2df * a1 = (v2df*) a[1];
    v2df * a2 = (v2df*) a[2];
    v2df * a3 = (v2df*) a[3];
    int nh = n>>1;
    int j;

    int k;
    v2df * cv;
    v2df * cvv;
    register v2df b0;
    register v2df b1;
    register v2df b2;
    register v2df b3;
    register v2df b4;
    register v2df b5;
    register v2df b6;
    register v2df b7;
    for(k=0;k<4;k+=2){
	cv = (v2df*) c[k];
	cvv = (v2df*) c[k+1];
	b0 = (v2df){b[k][0],b[k][0]};
	b1 = (v2df){b[k][1],b[k][1]};
	b2 = (v2df){b[k][2],b[k][2]};
	b3 = (v2df){b[k][3],b[k][3]};
	b4 = (v2df){b[k+1][0],b[k+1][0]};
	b5 = (v2df){b[k+1][1],b[k+1][1]};
	b6 = (v2df){b[k+1][2],b[k+1][2]};
	b7 = (v2df){b[k+1][3],b[k+1][3]};
	for(j=0;j<nh;j++){
	    register v2df aa0 = a0[j];
	    register v2df aa1 = a1[j];
	    register v2df aa2 = a2[j];
	    register v2df aa3 = a3[j];
	    register v2df x = aa0*b0;
	    x+= aa1*b1;
	    x+= aa2*b2;
	    x+= aa3*b3;
	    cv[j]-=x;
	    x = aa0*b4;
	    x+= aa1*b5;
	    x+= aa2*b6;
	    x+= aa3*b7;
	    cvv[j] -= x;
	}
    }
}

static void matmul_for_nk4(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int m0,m1,m2,m3;
    int s1, s2, s3;
    int ds = (n/32)*8;
    s1 = ds;
    s2 = ds*2;
    s3 = ds*3;
    m0 = ds;
    m1 = ds;
    m2 = ds;
    m3 = n-s3;
    //    fprintf(stderr,"n, s, m = %d %d %d %d %d %d %d %d\n",
    //    n,s1,s2,s3,m0,m1,m2,m3);

#pragma omp parallel
#pragma omp sections
  {
#pragma omp section
      matmul_for_nk4_7(n1,a,n2,b,n3,c,m0);
#pragma omp section
      matmul_for_nk4_7(n1,((double*)a)+s1 ,n2,b,n3,((double*)c)+s1,m1);
#pragma omp section
      matmul_for_nk4_7(n1,((double*)a)+s2 ,n2,b,n3,((double*)c)+s2,m2);
#pragma omp section
      matmul_for_nk4_7(n1,((double*)a)+s3 ,n2,b,n3,((double*)c)+s3,m3);
  }
}
static void matmul_for_nk8_9(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    matmul_for_nk4(n1, a, n2, b, n3,c, n);
    matmul_for_nk4(n1, (double(*)[])  (a[4]), n2,
		   (double(*)[])  &(b[0][4]), n3,c, n);
    matmul_for_nk4(n1, a, n2, (double(*)[])  b[4], n3,
		   (double(*)[])  c[4], n);
    matmul_for_nk4(n1, (double(*)[])  &(a[4]), n2,
		   (double(*)[])  &(b[4][4]), n3,
		   (double(*)[])  c[4], n);

}
#if 0
static void matmul_for_nk8(int n1, double a[][n1],
			   int n2, double b[][n2],
			   int n3, double c[][n3],
			   int n)
{
    int i;
    int nb=384;
    for (i=0;i<n;i+=nb){
	int iend = i+nb;
	if (iend > n) iend = n;
	matmul_for_nk8_worker(n1, (double(*)[]) (a[0]+i),
			      n2,  (double(*)[]) (b[0]+i),
			      n3, (double(*)[]) (c[0]+i),
			      iend-i);
    }
}
#endif

static void matmul_for_small_nk_local(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int kk,
			 int n)
{
    // simplest version
    int j,k,l;
    BEGIN_TSC;
    if (kk == 2){
	matmul_for_nk2(n1, a, n2, b, n3,c, n);
	END_TSC(t,16);
	return;
    }
    if (kk == 4){
	matmul_for_nk4(n1, a, n2, b, n3,c, n);
	END_TSC(t,13);
	return;
    }
    if (kk == 8){
	matmul_for_nk8(n1, a, n2, b, n3,c, n);
	END_TSC(t,12);
	return;
    }

    for(j=0;j<n;j++)
	for(k=0;k<m;k++)
	    for(l=0;l<kk;l++)
		c[k][j] -= a[l][j]*b[k][l];
}

static void mmmulandsub(int n, double a[n+1][n], int rshift, int m0, int m1,
		 int c0, int c1, int r0,int r1)
{
    int j,k,l;
#ifdef TIMETEST
    BEGIN_TSC;
#endif
#ifndef NOBLAS
    if ((m1-m0)<=8){
	// =4 is slighly faster than =8 on Ci7
	matmul_for_small_nk_local(n, (double(*)[])  &(a[m0-rshift][r0]),
			      n, (double(*)[])  &(a[c0-rshift][m0]),
			      n, (double(*)[])  &(a[c0-rshift][r0]),
			      c1-c0,m1-m0,r1-r0);
    }else{
	cblas_dgemm( CblasColMajor, CblasNoTrans, CblasNoTrans,
		     r1-r0, c1-c0, m1-m0, -1.0, &(a[m0-rshift][r0]), n,
		     &(a[c0-rshift][m0]), n, 1, &(a[c0-rshift][r0]), n );
    }
    // example:
    // r0, m0 = i+m,i
    // m0, c0 = i, i+m
    // r0, c0 = i+m, i+m
    //r1-r0 = n-i-m
    // c1-c0 = iend-i-m
    // m1-m0 = m
#else
    for(j=r0;j<r1;j++)
	for (k=c0;k<c1;k++)
	    for (l=m0; l<m1; l++)
		a[k-rshift][j] -= a[l-rshift][j]*a[k-rshift][l];
#endif
#ifdef TIMETEST
    END_TSC(t,4);
#endif
}

static int nswap;
static void column_decomposition(int n, double a[][n],  int m, int pv[], int i)
{
    // shift a so that partial array is okay
    int  j, k;
    int ip,ii;
    double ainv;
#ifdef TIMETEST
    BEGIN_TSC;
#endif
    for(ip=0;ip<m;ip++){
	ii=i+ip;
	int p = findpivot(n,a,ii);
	pv[ip]=p;
	swaprows(n,a,p,ii,0,m);
	nswap++;
	// normalize row ii
	ainv = 1.0/a[ip][ii];
	scalerow(n,a,ainv,ii,0,ip);
	scalerow(n,a,ainv,ii,ip+1,m);
	// subtract row ii from all lower rows
	vvmulandsub(n,  a, ip,ii, ip+1, m, ii+1, n);
    }
#ifdef TIMETEST
    END_TSC(t,5);
#endif
}

static void process_right_part(int n,
			double a[n+1][n],
			int m,
			int pv[],
			int i,
			int iend)
{
    int ii;
    // exchange rows
    for(ii=i;ii<i+m;ii++){
	swaprows(n,a,pv[ii-i],ii,m,iend-i);
    }

    // normalize rows
    for(ii=i;ii<i+m;ii++){
	scalerow(n,a,1.0/a[ii-i][ii] ,ii,m,iend-i);
    }
    // subtract rows (within i-i+m-1)
    for(ii=i;ii<i+m;ii++){
	vvmulandsub(n,  a, ii-i,ii,  m, iend-i, ii+1, i+m);
    }

    // subtract rows i-i+m-1 from all lower rows
    mmmulandsub(n, a, i, i,i+m, i+m, iend, i+m, n);
    //    fprintf(stderr,"process_r, end\n");
    //    usleep(1);
}
static void column_decomposition_recursive(int n,
				    double a[n+1][n],
				    int m,
				    int pv[],
				    int i)
{
    int  j, k;
    int ip,ii;
    double ainv;
    //    fprintf(stderr,"enter t column recursive %d %d\n", i, m);
    if (m <= 2){
	// perform non-recursive direct decomposition
	column_decomposition(n, a, m, pv,i);
    }else{
	// process the left half by recursion
	//	fprintf(stderr,"call column recursive %d %d\n", i, m);
	column_decomposition_recursive(n, a, m/2, pv,i);
	// process the right half
	//	fprintf(stderr,"call right part %d %d\n", i, m);
	process_right_part(n,a,m/2,pv,i,i+m);
	//	fprintf(stderr,"call right recursive %d %d\n", i, m);
	column_decomposition_recursive(n, a+m/2, m/2, pv+m/2,i+m/2);
	// process the swap of rows for the left half
	//	fprintf(stderr,"call swaprowse %d %d\n", i, m);
	for(ii=i+m/2;ii<i+m;ii++){
	    swaprows(n,a,pv[ii-i],ii,0,m/2);
	}
	//	fprintf(stderr,"call scalerows %d %d\n", i, m);
	// normalize rows
	for(ii=i+m/2;ii<i+m;ii++){
	    scalerow(n,a,1.0/a[ii-i][ii] ,ii,0,m/2);
	}
    }
}

void cm_column_decomposition_recursive(int n,
				    double a[n+1][n],
				    int m,
				    int pv[],
				    int i)
{
    column_decomposition_recursive( n, a, m, pv, i);
}
void cm_column_decomposition(int n,
				    double a[n+1][n],
				    int m,
				    int pv[],
				    int i)
{
    column_decomposition( n, a, m, pv, i);
}
void cm_process_right_part(int n,
			   double a[n+1][n],
			   int m,
			   int pv[],
			   int i,
			   int iend)
{
    process_right_part(n, a, m, pv, i,iend);
}


text:
/**
 *
 * @file example_strsmpl.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @brief Example for solving a system of linear equations using
 *        LU factorization and PLASMA_strsmpl routine
 *
 * @version 2.6.0
 * @author Bilel Hadri
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:45:21 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>

int check_solution(int, int , float *, int, float *, float *, int);

int IONE=1;
int ISEED[4] = {0,0,0,1};   /* initial seed for slarnv() */

int main ()
{

    int cores = 2;
    int N     = 10;
    int LDA   = 10;
    int NRHS  = 5;
    int LDB   = 10;
    int info;
    int info_solution;
    int i,j;
    int LDAxN = LDA*N;
    int LDBxNRHS = LDB*NRHS;

    float *A1 = (float *)malloc(LDA*N*(sizeof*A1));
    float *A2 = (float *)malloc(LDA*N*(sizeof*A2));
    float *B1 = (float *)malloc(LDB*NRHS*(sizeof*B1));
    float *B2 = (float *)malloc(LDB*NRHS*(sizeof*B2));
    PLASMA_desc *L;
    int *IPIV;

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)||(!B1)||(!B2)){
        printf("Out of Memory \n ");
        return EXIT_SUCCESS;
    }

    /*Plasma Initialize*/
    PLASMA_Init(cores);
    printf("-- PLASMA is initialized to run on %d cores. \n",cores);

    /* Initialize A1 and A2 Matrix */
    LAPACKE_slarnv_work(IONE, ISEED, LDAxN, A1);
    for ( i = 0; i < N; i++)
        for (  j = 0; j < N; j++)
            A2[LDA*j+i] = A1[LDA*j+i];

    /* Initialize B1 and B2 */
    LAPACKE_slarnv_work(IONE, ISEED, LDBxNRHS, B1);
    for ( i = 0; i < N; i++)
        for ( j = 0; j < NRHS; j++)
            B2[LDB*j+i] = B1[LDB*j+i];

    /* Allocate L and IPIV */
    info = PLASMA_Alloc_Workspace_sgetrf_incpiv(N, N, &L, &IPIV);

    /* LU factorization of the matrix A */
    info = PLASMA_sgetrf_incpiv(N, N, A2, LDA, L, IPIV);

    /* Solve the problem */
    info = PLASMA_strsmpl(N, NRHS, A2, LDA, L, IPIV, B2, LDB);
    info = PLASMA_strsm(PlasmaLeft, PlasmaUpper, PlasmaNoTrans, PlasmaNonUnit,
                        N, NRHS, (float)1.0, A2, LDA, B2, LDB);

    /* Check the solution */
    info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB);

    if ((info_solution != 0)|(info != 0))
       printf("-- Error in SGETRS example ! \n");
    else
       printf("-- Run of SGETRS example successful ! \n");

    free(A1); free(A2); free(B1); free(B2); free(IPIV); free(L);

    PLASMA_Finalize();

    return EXIT_SUCCESS;
}

/*------------------------------------------------------------------------
 *  Check the accuracy of the solution of the linear system
 */

int check_solution(int N, int NRHS, float *A1, int LDA, float *B1, float *B2, int LDB)
{
    int info_solution;
    float Rnorm, Anorm, Xnorm, Bnorm;
    float alpha, beta;
    float *work = (float *)malloc(N*sizeof(float));
    float eps;

    eps = LAPACKE_slamch_work('e');

    alpha = 1.0;
    beta  = -1.0;

    Xnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, NRHS, B2, LDB, work);
    Anorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, A1, LDA, work);
    Bnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, NRHS, B1, LDB, work);

    cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, NRHS, N, (alpha), A1, LDA, B2, LDB, (beta), B1, LDB);
    Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, NRHS, B1, LDB, work);

    printf("============\n");
    printf("Checking the Residual of the solution \n");
    printf("-- ||Ax-B||_oo/((||A||_oo||x||_oo+||B||_oo).N.eps) = %e \n",Rnorm/((Anorm*Xnorm+Bnorm)*N*eps));

    if ( isnan(Rnorm/((Anorm*Xnorm+Bnorm)*N*eps)) || (Rnorm/((Anorm*Xnorm+Bnorm)*N*eps) > 10.0) ){
        printf("-- The solution is suspicious ! \n");
        info_solution = 1;
    }
    else{
        printf("-- The solution is CORRECT ! \n");
        info_solution = 0;
    }

    free(work);

    return info_solution;
}

ext: c


text:
#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>

int main(void)
{
  double x = 5.0;
  double y = gsl_sf_bessel_J0(x);
  printf("J0(%g) = %.18e\n", x, y);
  return 0;
}

max_stars_repo_name: jbelmont/algorithms-workshop

max_issues_repo_name: jbelmont/algorithms-workshop

max_forks_repo_name: jbelmont/algorithms-workshop


text:
#pragma once

#include <emerald/util/foundation.h>

#define gsl_FEATURE_MAKE_SPAN_TO_STD 20
#include <gsl/gsl-lite.hpp>

#include <array>

namespace emerald::shallow_weaver {

using namespace emerald::util;

template <typename T>
using span = gsl::span<T>;

enum class Time_integration { FIRST_ORDER, RUNGE_KUTTA_2, RUNGE_KUTTA_4 };

using int2 = std::array<int, 2>;
using int4 = std::array<int, 4>;
using float2 = std::array<float, 2>;
using float4 = std::array<float, 4>;

}  // namespace emerald::shallow_weaver

ext: h

max_stars_repo_name: blackencino/emerald

max_issues_repo_name: blackencino/emerald

max_forks_repo_name: blackencino/emerald


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * 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 2 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, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:43:58 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -real2hc 16 */

/*
 * This function contains 58 FP additions, 12 FP multiplications,
 * (or, 54 additions, 8 multiplications, 4 fused multiply/add),
 * 30 stack variables, and 32 memory accesses
 */
static const fftw_real K707106781 = FFTW_KONST(+0.707106781186547524400844362104849039284835938);
static const fftw_real K923879532 = FFTW_KONST(+0.923879532511286756128183189396788286822416626);
static const fftw_real K382683432 = FFTW_KONST(+0.382683432365089771728459984030398866761344562);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftw_real2hc_16(const fftw_real *input, fftw_real *real_output, fftw_real *imag_output, int istride, int real_ostride, int imag_ostride)
{
     fftw_real tmp3;
     fftw_real tmp6;
     fftw_real tmp7;
     fftw_real tmp35;
     fftw_real tmp18;
     fftw_real tmp33;
     fftw_real tmp40;
     fftw_real tmp48;
     fftw_real tmp56;
     fftw_real tmp10;
     fftw_real tmp13;
     fftw_real tmp14;
     fftw_real tmp36;
     fftw_real tmp17;
     fftw_real tmp26;
     fftw_real tmp41;
     fftw_real tmp51;
     fftw_real tmp57;
     fftw_real tmp16;
     fftw_real tmp15;
     fftw_real tmp43;
     fftw_real tmp44;
     ASSERT_ALIGNED_DOUBLE;
     {
	  fftw_real tmp1;
	  fftw_real tmp2;
	  fftw_real tmp4;
	  fftw_real tmp5;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp1 = input[0];
	  tmp2 = input[8 * istride];
	  tmp3 = tmp1 + tmp2;
	  tmp4 = input[4 * istride];
	  tmp5 = input[12 * istride];
	  tmp6 = tmp4 + tmp5;
	  tmp7 = tmp3 + tmp6;
	  tmp35 = tmp1 - tmp2;
	  tmp18 = tmp4 - tmp5;
     }
     {
	  fftw_real tmp29;
	  fftw_real tmp46;
	  fftw_real tmp32;
	  fftw_real tmp47;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp27;
	       fftw_real tmp28;
	       fftw_real tmp30;
	       fftw_real tmp31;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp27 = input[istride];
	       tmp28 = input[9 * istride];
	       tmp29 = tmp27 - tmp28;
	       tmp46 = tmp27 + tmp28;
	       tmp30 = input[5 * istride];
	       tmp31 = input[13 * istride];
	       tmp32 = tmp30 - tmp31;
	       tmp47 = tmp30 + tmp31;
	  }
	  tmp33 = (K382683432 * tmp29) + (K923879532 * tmp32);
	  tmp40 = (K923879532 * tmp29) - (K382683432 * tmp32);
	  tmp48 = tmp46 - tmp47;
	  tmp56 = tmp46 + tmp47;
     }
     {
	  fftw_real tmp8;
	  fftw_real tmp9;
	  fftw_real tmp11;
	  fftw_real tmp12;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp8 = input[2 * istride];
	  tmp9 = input[10 * istride];
	  tmp10 = tmp8 + tmp9;
	  tmp16 = tmp8 - tmp9;
	  tmp11 = input[14 * istride];
	  tmp12 = input[6 * istride];
	  tmp13 = tmp11 + tmp12;
	  tmp15 = tmp11 - tmp12;
     }
     tmp14 = tmp10 + tmp13;
     tmp36 = K707106781 * (tmp16 + tmp15);
     tmp17 = K707106781 * (tmp15 - tmp16);
     {
	  fftw_real tmp22;
	  fftw_real tmp49;
	  fftw_real tmp25;
	  fftw_real tmp50;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp20;
	       fftw_real tmp21;
	       fftw_real tmp23;
	       fftw_real tmp24;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp20 = input[15 * istride];
	       tmp21 = input[7 * istride];
	       tmp22 = tmp20 - tmp21;
	       tmp49 = tmp20 + tmp21;
	       tmp23 = input[3 * istride];
	       tmp24 = input[11 * istride];
	       tmp25 = tmp23 - tmp24;
	       tmp50 = tmp23 + tmp24;
	  }
	  tmp26 = (K382683432 * tmp22) - (K923879532 * tmp25);
	  tmp41 = (K923879532 * tmp22) + (K382683432 * tmp25);
	  tmp51 = tmp49 - tmp50;
	  tmp57 = tmp49 + tmp50;
     }
     {
	  fftw_real tmp55;
	  fftw_real tmp58;
	  fftw_real tmp53;
	  fftw_real tmp54;
	  ASSERT_ALIGNED_DOUBLE;
	  real_output[4 * real_ostride] = tmp7 - tmp14;
	  tmp55 = tmp7 + tmp14;
	  tmp58 = tmp56 + tmp57;
	  real_output[8 * real_ostride] = tmp55 - tmp58;
	  real_output[0] = tmp55 + tmp58;
	  imag_output[4 * imag_ostride] = tmp57 - tmp56;
	  tmp53 = tmp13 - tmp10;
	  tmp54 = K707106781 * (tmp51 - tmp48);
	  imag_output[2 * imag_ostride] = tmp53 + tmp54;
	  imag_output[6 * imag_ostride] = tmp54 - tmp53;
     }
     {
	  fftw_real tmp45;
	  fftw_real tmp52;
	  fftw_real tmp39;
	  fftw_real tmp42;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp45 = tmp3 - tmp6;
	  tmp52 = K707106781 * (tmp48 + tmp51);
	  real_output[6 * real_ostride] = tmp45 - tmp52;
	  real_output[2 * real_ostride] = tmp45 + tmp52;
	  tmp39 = tmp35 + tmp36;
	  tmp42 = tmp40 + tmp41;
	  real_output[7 * real_ostride] = tmp39 - tmp42;
	  real_output[real_ostride] = tmp39 + tmp42;
     }
     tmp43 = tmp18 + tmp17;
     tmp44 = tmp41 - tmp40;
     imag_output[3 * imag_ostride] = tmp43 + tmp44;
     imag_output[5 * imag_ostride] = tmp44 - tmp43;
     {
	  fftw_real tmp19;
	  fftw_real tmp34;
	  fftw_real tmp37;
	  fftw_real tmp38;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp19 = tmp17 - tmp18;
	  tmp34 = tmp26 - tmp33;
	  imag_output[imag_ostride] = tmp19 + tmp34;
	  imag_output[7 * imag_ostride] = tmp34 - tmp19;
	  tmp37 = tmp35 - tmp36;
	  tmp38 = tmp33 + tmp26;
	  real_output[5 * real_ostride] = tmp37 - tmp38;
	  real_output[3 * real_ostride] = tmp37 + tmp38;
     }
}

fftw_codelet_desc fftw_real2hc_16_desc =
{
     "fftw_real2hc_16",
     (void (*)()) fftw_real2hc_16,
     16,
     FFTW_FORWARD,
     FFTW_REAL2HC,
     354,
     0,
     (const int *) 0,
};

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt


text:
/* specfunc/legendre_Qn.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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 2 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, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_bessel.h"
#include "gsl_sf_elementary.h"
#include "gsl_sf_exp.h"
#include "gsl_sf_pow_int.h"
#include "gsl_sf_legendre.h"

#include "error.h"

/* Evaluate f_{ell+1}/f_ell
 * f_ell := Q^{b}_{a+ell}(x)
 * x > 1
 */
static
int
legendreQ_CF1_xgt1(int ell, double a, double b, double x, double * result)
{
  const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  const int maxiter = 5000;
  int n = 1;
  double Anm2 = 1.0;
  double Bnm2 = 0.0;
  double Anm1 = 0.0;
  double Bnm1 = 1.0;
  double a1 = ell + 1.0 + a + b;
  double b1 = (2.0*(ell+1.0+a) + 1.0) * x;
  double An = b1*Anm1 + a1*Anm2;
  double Bn = b1*Bnm1 + a1*Bnm2;
  double an, bn;
  double fn = An/Bn;

  while(n < maxiter) {
    double old_fn;
    double del;
    double lna;
    n++;
    Anm2 = Anm1;
    Bnm2 = Bnm1;
    Anm1 = An;
    Bnm1 = Bn;
    lna = ell + n + a;
    an = b*b - lna*lna;
    bn = (2.0*lna + 1.0) * x;
    An = bn*Anm1 + an*Anm2;
    Bn = bn*Bnm1 + an*Bnm2;

    if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
      An /= RECUR_BIG;
      Bn /= RECUR_BIG;
      Anm1 /= RECUR_BIG;
      Bnm1 /= RECUR_BIG;
      Anm2 /= RECUR_BIG;
      Bnm2 /= RECUR_BIG;
    }

    old_fn = fn;
    fn = An/Bn;
    del = old_fn/fn;

    if(fabs(del - 1.0) < 4.0*GSL_DBL_EPSILON) break;
  }

  *result = fn;

  if(n == maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}

/* Uniform asymptotic for Q_l(x).
 * Assumes x > -1.0 and x != 1.0.
 * Discards second order and higher terms.
 */
static
int
legendre_Ql_asymp_unif(const double ell, const double x, gsl_sf_result * result)
{
  if(x < 1.0) {
    double u   = ell + 0.5;
    double th  = acos(x);
    gsl_sf_result Y0, Y1;
    int stat_Y0, stat_Y1;
    int stat_m;
    double pre;
    double B00;
    double sum;

    /* B00 = 1/8 (1 - th cot(th) / th^2
     * pre = sqrt(th/sin(th))
     */
    if(th < GSL_ROOT4_DBL_EPSILON) {
      B00 = (1.0 + th*th/15.0)/24.0;
      pre = 1.0 + th*th/12.0;
    }
    else {
      double sin_th = sqrt(1.0 - x*x);
      double cot_th = x / sin_th;
      B00 = 1.0/8.0 * (1.0 - th * cot_th) / (th*th);
      pre = sqrt(th/sin_th);
    }

    stat_Y0 = gsl_sf_bessel_Y0_e(u*th, &Y0);
    stat_Y1 = gsl_sf_bessel_Y1_e(u*th, &Y1);

    sum = -0.5*M_PI * (Y0.val + th/u * Y1.val * B00);

    stat_m = gsl_sf_multiply_e(pre, sum, result);
    result->err += 0.5*M_PI * fabs(pre) * (Y0.err + fabs(th/u*B00)*Y1.err);
    result->err += GSL_DBL_EPSILON * fabs(result->val);

    return GSL_ERROR_SELECT_3(stat_m, stat_Y0, stat_Y1);
  }
  else {
    double u   = ell + 0.5;
    double xi  = acosh(x);
    gsl_sf_result K0_scaled, K1_scaled;
    int stat_K0, stat_K1;
    int stat_e;
    double pre;
    double B00;
    double sum;

    /* B00 = -1/8 (1 - xi coth(xi) / xi^2
     * pre = sqrt(xi/sinh(xi))
     */
    if(xi < GSL_ROOT4_DBL_EPSILON) {
      B00 = (1.0-xi*xi/15.0)/24.0;
      pre = 1.0 - xi*xi/12.0;
    }
    else {
      double sinh_xi = sqrt(x*x - 1.0);
      double coth_xi = x / sinh_xi;
      B00 = -1.0/8.0 * (1.0 - xi * coth_xi) / (xi*xi);
      pre = sqrt(xi/sinh_xi);
    }

    stat_K0 = gsl_sf_bessel_K0_scaled_e(u*xi, &K0_scaled);
    stat_K1 = gsl_sf_bessel_K1_scaled_e(u*xi, &K1_scaled);

    sum = K0_scaled.val - xi/u * K1_scaled.val * B00;

    stat_e = gsl_sf_exp_mult_e(-u*xi, pre * sum, result);
    result->err  = GSL_DBL_EPSILON * fabs(result->val) * fabs(u*xi);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_ERROR_SELECT_3(stat_e, stat_K0, stat_K1);
  }
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_legendre_Q0_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= -1.0 || x == 1.0) {
    DOMAIN_ERROR(result);
  }
  else if(x*x < GSL_ROOT6_DBL_EPSILON) { /* |x| <~ 0.05 */
    const double c3 = 1.0/3.0;
    const double c5 = 1.0/5.0;
    const double c7 = 1.0/7.0;
    const double c9 = 1.0/9.0;
    const double c11 = 1.0/11.0;
    const double y = x * x;
    const double series = 1.0 + y*(c3 + y*(c5 + y*(c7 + y*(c9 + y*c11))));
    result->val = x * series;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(x);
    return GSL_SUCCESS;
  }
  else if(x < 1.0) {
    result->val = 0.5 * log((1.0+x)/(1.0-x));
    result->err  = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < 10.0) {
    result->val = 0.5 * log((x+1.0)/(x-1.0));
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x*GSL_DBL_MIN < 2.0) {
    const double y = 1.0/(x*x);
    const double c1 = 1.0/3.0;
    const double c2 = 1.0/5.0;
    const double c3 = 1.0/7.0;
    const double c4 = 1.0/9.0;
    const double c5 = 1.0/11.0;
    const double c6 = 1.0/13.0;
    const double c7 = 1.0/15.0;
    result->val = (1.0/x) * (1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*(c5 + y*(c6 + y*c7)))))));
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    UNDERFLOW_ERROR(result);
  }
}

int
gsl_sf_legendre_Q1_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= -1.0 || x == 1.0) {
    DOMAIN_ERROR(result);
  }
  else if(x*x < GSL_ROOT6_DBL_EPSILON) { /* |x| <~ 0.05 */
    const double c3 = 1.0/3.0;
    const double c5 = 1.0/5.0;
    const double c7 = 1.0/7.0;
    const double c9 = 1.0/9.0;
    const double c11 = 1.0/11.0;
    const double y = x * x;
    const double series = 1.0 + y*(c3 + y*(c5 + y*(c7 + y*(c9 + y*c11))));
    result->val = x * x * series - 1.0;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < 1.0){
    result->val = 0.5 * x * (log((1.0+x)/(1.0-x))) - 1.0;
    result->err  = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < 6.0) {
    result->val = 0.5 * x * log((x+1.0)/(x-1.0)) - 1.0;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x*GSL_SQRT_DBL_MIN < 0.99/M_SQRT3) {
    const double y = 1/(x*x);
    const double c1 = 3.0/5.0;
    const double c2 = 3.0/7.0;
    const double c3 = 3.0/9.0;
    const double c4 = 3.0/11.0;
    const double c5 = 3.0/13.0;
    const double c6 = 3.0/15.0;
    const double c7 = 3.0/17.0;
    const double c8 = 3.0/19.0;
    const double sum = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*(c5 + y*(c6 + y*(c7 + y*c8)))))));
    result->val = sum / (3.0*x*x);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    UNDERFLOW_ERROR(result);
  }
}

int
gsl_sf_legendre_Ql_e(const int l, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= -1.0 || x == 1.0 || l < 0) {
    DOMAIN_ERROR(result);
  }
  else if(l == 0) {
    return gsl_sf_legendre_Q0_e(x, result);
  }
  else if(l == 1) {
    return gsl_sf_legendre_Q1_e(x, result);
  }
  else if(l > 100000) {
    return legendre_Ql_asymp_unif(l, x, result);
  }
  else if(x < 1.0){
    /* Forward recurrence.
     */
    gsl_sf_result Q0, Q1;
    int stat_Q0 = gsl_sf_legendre_Q0_e(x, &Q0);
    int stat_Q1 = gsl_sf_legendre_Q1_e(x, &Q1);
    double Qellm1 = Q0.val;
    double Qell   = Q1.val;
    double Qellp1;
    int ell;
    for(ell=1; ell<l; ell++) {
      Qellp1 = (x*(2.0*ell + 1.0) * Qell - ell * Qellm1) / (ell + 1.0);
      Qellm1 = Qell;
      Qell   = Qellp1;
    }
    result->val = Qell;
    result->err = GSL_DBL_EPSILON * l * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_Q0, stat_Q1);
  }
  else {
    /* x > 1.0 */

    double rat;
    int stat_CF1  = legendreQ_CF1_xgt1(l, 0.0, 0.0, x, &rat);
    int stat_Q;
    double Qellp1 = rat * GSL_SQRT_DBL_MIN;
    double Qell   = GSL_SQRT_DBL_MIN;
    double Qellm1;
    int ell;
    for(ell=l; ell>0; ell--) {
      Qellm1 = (x * (2.0*ell + 1.0) * Qell - (ell+1.0) * Qellp1) / ell;
      Qellp1 = Qell;
      Qell   = Qellm1;
    }

    if(fabs(Qell) > fabs(Qellp1)) {
      gsl_sf_result Q0;
      stat_Q = gsl_sf_legendre_Q0_e(x, &Q0);
      result->val = GSL_SQRT_DBL_MIN * Q0.val / Qell;
      result->err = l * GSL_DBL_EPSILON * fabs(result->val);
    }
    else {
      gsl_sf_result Q1;
      stat_Q = gsl_sf_legendre_Q1_e(x, &Q1);
      result->val = GSL_SQRT_DBL_MIN * Q1.val / Qellp1;
      result->err = l * GSL_DBL_EPSILON * fabs(result->val);
    }

    return GSL_ERROR_SELECT_2(stat_Q, stat_CF1);
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_legendre_Q0(const double x)
{
  EVAL_RESULT(gsl_sf_legendre_Q0_e(x, &result));
}

double gsl_sf_legendre_Q1(const double x)
{
  EVAL_RESULT(gsl_sf_legendre_Q1_e(x, &result));
}

double gsl_sf_legendre_Ql(const int l, const double x)
{
  EVAL_RESULT(gsl_sf_legendre_Ql_e(l, x, &result));
}


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * 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 2 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, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:44:26 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -hc2real 15 */

/*
 * This function contains 64 FP additions, 31 FP multiplications,
 * (or, 47 additions, 14 multiplications, 17 fused multiply/add),
 * 36 stack variables, and 30 memory accesses
 */
static const fftw_real K1_118033988 = FFTW_KONST(+1.118033988749894848204586834365638117720309180);
static const fftw_real K1_902113032 = FFTW_KONST(+1.902113032590307144232878666758764286811397268);
static const fftw_real K1_175570504 = FFTW_KONST(+1.175570504584946258337411909278145537195304875);
static const fftw_real K500000000 = FFTW_KONST(+0.500000000000000000000000000000000000000000000);
static const fftw_real K866025403 = FFTW_KONST(+0.866025403784438646763723170752936183471402627);
static const fftw_real K2_000000000 = FFTW_KONST(+2.000000000000000000000000000000000000000000000);
static const fftw_real K1_732050807 = FFTW_KONST(+1.732050807568877293527446341505872366942805254);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftw_hc2real_15(const fftw_real *real_input, const fftw_real *imag_input, fftw_real *output, int real_istride, int imag_istride, int ostride)
{
     fftw_real tmp3;
     fftw_real tmp30;
     fftw_real tmp18;
     fftw_real tmp37;
     fftw_real tmp61;
     fftw_real tmp62;
     fftw_real tmp45;
     fftw_real tmp40;
     fftw_real tmp23;
     fftw_real tmp31;
     fftw_real tmp42;
     fftw_real tmp28;
     fftw_real tmp32;
     fftw_real tmp8;
     fftw_real tmp13;
     fftw_real tmp14;
     ASSERT_ALIGNED_DOUBLE;
     {
	  fftw_real tmp17;
	  fftw_real tmp1;
	  fftw_real tmp2;
	  fftw_real tmp15;
	  fftw_real tmp16;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp16 = imag_input[5 * imag_istride];
	  tmp17 = K1_732050807 * tmp16;
	  tmp1 = real_input[0];
	  tmp2 = real_input[5 * real_istride];
	  tmp15 = tmp1 - tmp2;
	  tmp3 = tmp1 + (K2_000000000 * tmp2);
	  tmp30 = tmp15 - tmp17;
	  tmp18 = tmp15 + tmp17;
     }
     {
	  fftw_real tmp4;
	  fftw_real tmp39;
	  fftw_real tmp5;
	  fftw_real tmp6;
	  fftw_real tmp7;
	  fftw_real tmp22;
	  fftw_real tmp38;
	  fftw_real tmp9;
	  fftw_real tmp44;
	  fftw_real tmp10;
	  fftw_real tmp11;
	  fftw_real tmp12;
	  fftw_real tmp27;
	  fftw_real tmp43;
	  fftw_real tmp19;
	  fftw_real tmp24;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp20;
	       fftw_real tmp21;
	       fftw_real tmp25;
	       fftw_real tmp26;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp4 = real_input[3 * real_istride];
	       tmp39 = imag_input[3 * imag_istride];
	       tmp5 = real_input[7 * real_istride];
	       tmp6 = real_input[2 * real_istride];
	       tmp7 = tmp5 + tmp6;
	       tmp20 = imag_input[7 * imag_istride];
	       tmp21 = imag_input[2 * imag_istride];
	       tmp22 = K866025403 * (tmp20 - tmp21);
	       tmp38 = tmp20 + tmp21;
	       tmp9 = real_input[6 * real_istride];
	       tmp44 = imag_input[6 * imag_istride];
	       tmp10 = real_input[4 * real_istride];
	       tmp11 = real_input[real_istride];
	       tmp12 = tmp10 + tmp11;
	       tmp25 = imag_input[4 * imag_istride];
	       tmp26 = imag_input[imag_istride];
	       tmp27 = K866025403 * (tmp25 + tmp26);
	       tmp43 = tmp25 - tmp26;
	  }
	  tmp37 = K866025403 * (tmp5 - tmp6);
	  tmp61 = tmp39 - tmp38;
	  tmp62 = tmp44 - tmp43;
	  tmp45 = (K500000000 * tmp43) + tmp44;
	  tmp40 = (K500000000 * tmp38) + tmp39;
	  tmp19 = tmp4 - (K500000000 * tmp7);
	  tmp23 = tmp19 - tmp22;
	  tmp31 = tmp19 + tmp22;
	  tmp42 = K866025403 * (tmp10 - tmp11);
	  tmp24 = tmp9 - (K500000000 * tmp12);
	  tmp28 = tmp24 - tmp27;
	  tmp32 = tmp24 + tmp27;
	  tmp8 = tmp4 + tmp7;
	  tmp13 = tmp9 + tmp12;
	  tmp14 = tmp8 + tmp13;
     }
     output[0] = tmp3 + (K2_000000000 * tmp14);
     {
	  fftw_real tmp63;
	  fftw_real tmp65;
	  fftw_real tmp60;
	  fftw_real tmp64;
	  fftw_real tmp58;
	  fftw_real tmp59;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp63 = (K1_175570504 * tmp61) - (K1_902113032 * tmp62);
	  tmp65 = (K1_902113032 * tmp61) + (K1_175570504 * tmp62);
	  tmp58 = tmp3 - (K500000000 * tmp14);
	  tmp59 = K1_118033988 * (tmp8 - tmp13);
	  tmp60 = tmp58 - tmp59;
	  tmp64 = tmp59 + tmp58;
	  output[12 * ostride] = tmp60 - tmp63;
	  output[3 * ostride] = tmp60 + tmp63;
	  output[6 * ostride] = tmp64 - tmp65;
	  output[9 * ostride] = tmp64 + tmp65;
     }
     {
	  fftw_real tmp51;
	  fftw_real tmp29;
	  fftw_real tmp50;
	  fftw_real tmp55;
	  fftw_real tmp57;
	  fftw_real tmp53;
	  fftw_real tmp54;
	  fftw_real tmp56;
	  fftw_real tmp52;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp51 = K1_118033988 * (tmp23 - tmp28);
	  tmp29 = tmp23 + tmp28;
	  tmp50 = tmp18 - (K500000000 * tmp29);
	  tmp53 = tmp40 - tmp37;
	  tmp54 = tmp45 - tmp42;
	  tmp55 = (K1_175570504 * tmp53) - (K1_902113032 * tmp54);
	  tmp57 = (K1_902113032 * tmp53) + (K1_175570504 * tmp54);
	  output[5 * ostride] = tmp18 + (K2_000000000 * tmp29);
	  tmp56 = tmp51 + tmp50;
	  output[11 * ostride] = tmp56 - tmp57;
	  output[14 * ostride] = tmp56 + tmp57;
	  tmp52 = tmp50 - tmp51;
	  output[2 * ostride] = tmp52 - tmp55;
	  output[8 * ostride] = tmp52 + tmp55;
     }
     {
	  fftw_real tmp35;
	  fftw_real tmp33;
	  fftw_real tmp34;
	  fftw_real tmp47;
	  fftw_real tmp49;
	  fftw_real tmp41;
	  fftw_real tmp46;
	  fftw_real tmp48;
	  fftw_real tmp36;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp35 = K1_118033988 * (tmp31 - tmp32);
	  tmp33 = tmp31 + tmp32;
	  tmp34 = tmp30 - (K500000000 * tmp33);
	  tmp41 = tmp37 + tmp40;
	  tmp46 = tmp42 + tmp45;
	  tmp47 = (K1_175570504 * tmp41) - (K1_902113032 * tmp46);
	  tmp49 = (K1_902113032 * tmp41) + (K1_175570504 * tmp46);
	  output[10 * ostride] = tmp30 + (K2_000000000 * tmp33);
	  tmp48 = tmp35 + tmp34;
	  output[ostride] = tmp48 - tmp49;
	  output[4 * ostride] = tmp48 + tmp49;
	  tmp36 = tmp34 - tmp35;
	  output[7 * ostride] = tmp36 - tmp47;
	  output[13 * ostride] = tmp36 + tmp47;
     }
}

fftw_codelet_desc fftw_hc2real_15_desc =
{
     "fftw_hc2real_15",
     (void (*)()) fftw_hc2real_15,
     15,
     FFTW_BACKWARD,
     FFTW_HC2REAL,
     345,
     0,
     (const int *) 0,
};


text:
/* dht/dht.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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 2 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, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_bessel.h>
#include "gsl_dht.h"

gsl_dht *
gsl_dht_alloc (size_t size)
{
  gsl_dht * t;

  if(size == 0) {
    GSL_ERROR_VAL("size == 0", GSL_EDOM, 0);
  }

  t = (gsl_dht *)malloc(sizeof(gsl_dht));

  if(t == 0) {
    GSL_ERROR_VAL("out of memory", GSL_ENOMEM, 0);
  }

  t->size = size;

  t->xmax = -1.0; /* Make it clear that this needs to be calculated. */
  t->nu   = -1.0;

  t->j = (double *)malloc((size+2)*sizeof(double));

  if(t->j == 0) {
    free(t);
    GSL_ERROR_VAL("could not allocate memory for j", GSL_ENOMEM, 0);
  }

  t->Jjj = (double *)malloc(size*(size+1)/2 * sizeof(double));

  if(t->Jjj == 0) {
    free(t->j);
    free(t);
    GSL_ERROR_VAL("could not allocate memory for Jjj", GSL_ENOMEM, 0);
  }

  t->J2 = (double *)malloc((size+1)*sizeof(double));

  if(t->J2 == 0) {
    free(t->Jjj);
    free(t->j);
    free(t);
    GSL_ERROR_VAL("could not allocate memory for J2", GSL_ENOMEM, 0);
  }

  return t;
}

/* Handle internal calculation of Bessel zeros. */
static int
dht_bessel_zeros(gsl_dht * t)
{
  unsigned int s;
  gsl_sf_result z;
  int stat_z = 0;
  t->j[0] = 0.0;
  for(s=1; s < t->size + 2; s++) {
    stat_z += gsl_sf_bessel_zero_Jnu_e(t->nu, s, &z);
    t->j[s] = z.val;
  }
  if(stat_z != 0) {
    GSL_ERROR("could not compute bessel zeroes", GSL_EFAILED);
  }
  else {
    return GSL_SUCCESS;
  }
}

gsl_dht *
gsl_dht_new (size_t size, double nu, double xmax)
{
  int status;

  gsl_dht * dht = gsl_dht_alloc (size);

  if (dht == 0)
    return 0;

  status = gsl_dht_init(dht, nu, xmax);

  if (status)
    return 0;

  return dht;
}

int
gsl_dht_init(gsl_dht * t, double nu, double xmax)
{
  if(xmax <= 0.0) {
    GSL_ERROR ("xmax is not positive", GSL_EDOM);
  } else if(nu < 0.0) {
    GSL_ERROR ("nu is negative", GSL_EDOM);
  }
  else {
    size_t n, m;
    int stat_bz = GSL_SUCCESS;
    int stat_J  = 0;
    double jN;

    if(nu != t->nu) {
      /* Recalculate Bessel zeros if necessary. */
      t->nu = nu;
      stat_bz = dht_bessel_zeros(t);
    }

    jN = t->j[t->size+1];

    t->xmax = xmax;
    t->kmax = jN / xmax;

    t->J2[0] = 0.0;
    for(m=1; m<t->size+1; m++) {
      gsl_sf_result J;
      stat_J += gsl_sf_bessel_Jnu_e(nu + 1.0, t->j[m], &J);
      t->J2[m] = J.val * J.val;
    }

    /* J_nu(j[n] j[m] / j[N]) = Jjj[n(n-1)/2 + m - 1], 1 <= n,m <= size
     */
    for(n=1; n<t->size+1; n++) {
      for(m=1; m<=n; m++) {
        double arg = t->j[n] * t->j[m] / jN;
        gsl_sf_result J;
        stat_J += gsl_sf_bessel_Jnu_e(nu, arg, &J);
        t->Jjj[n*(n-1)/2 + m - 1] = J.val;
      }
    }

    if(stat_J != 0) {
      GSL_ERROR("error computing bessel function", GSL_EFAILED);
    }
    else {
      return stat_bz;
    }
  }
}

double gsl_dht_x_sample(const gsl_dht * t, int n)
{
  return t->j[n+1]/t->j[t->size+1] * t->xmax;
}

double gsl_dht_k_sample(const gsl_dht * t, int n)
{
  return t->j[n+1] / t->xmax;
}

void gsl_dht_free(gsl_dht * t)
{
  free(t->J2);
  free(t->Jjj);
  free(t->j);
  free(t);
}

int
gsl_dht_apply(const gsl_dht * t, double * f_in, double * f_out)
{
  const double jN = t->j[t->size + 1];
  const double r  = t->xmax / jN;
  size_t m;
  size_t i;

  for(m=0; m<t->size; m++) {
    double sum = 0.0;
    double Y;
    for(i=0; i<t->size; i++) {
      /* Need to find max and min so that we
       * address the symmetric Jjj matrix properly.
       * FIXME: we can presumably optimize this
       * by just running over the elements of Jjj
       * in a deterministic manner.
       */
      size_t m_local;
      size_t n_local;
      if(i < m) {
        m_local = i;
	n_local = m;
      }
      else {
        m_local = m;
	n_local = i;
      }
      Y = t->Jjj[n_local*(n_local+1)/2 + m_local] / t->J2[i+1];
      sum += Y * f_in[i];
    }
    f_out[m] = sum * 2.0 * r*r;
  }

  return GSL_SUCCESS;
}