text:
//
// lu2lib.c
//
// J. Makino
//    Time-stamp: <2019-05-06 22:08:30 makino>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <timerlib.h>
#ifndef NOBLAS
#ifdef MKL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#endif
#ifdef USEGDR
#include "gdrdgemm.h"
#endif

#define FTYPE double
#include <emmintrin.h>
typedef double v2df __attribute__((vector_size(16)));
typedef union {v2df v; double s[2];}v2u;

#ifndef USEGDR
void  gdrsetboardid(int boardid)
{}
#endif

void matmul2_host(int n,
		  FTYPE a[n][n],
		  FTYPE b[n][n],
		  FTYPE c[n][n])

{
    int i, j, k;
    for(i=0;i<n;i++){
	for(j=0;j<n;j++){
	    c[i][j]=0.0e0;
	    for(k=0;k<n;k++) c[i][j]+= a[i][k]*b[k][j];
	}
    }
}

// simplest version
void matmul_for_small_nk_0(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int kk,
			 int n)
{
    // simplest version
    int i,j,k;
    for(j=0;j<n;j++)
	for(i=0;i<m;i++)
	    for(k=0;k<kk;k++)
		c[i][j] -= a[i][k]*b[k][j];
}

// make copy of B
void matmul_for_small_nk_1(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int kk,
			 int n)
{

    int i,j,k;
    double bcopy[n][kk];
    for(k=0;k<kk;k++)
	for(j=0;j<n;j++)
		bcopy[j][k] = b[k][j];
    for(i=0;i<m;i++){
	for(j=0;j<n;j++){
	    register double tmp=0.0;
	    for(k=0;k<kk;k++){
		tmp += a[i][k]*bcopy[j][k];
	    }
	    c[i][j] -= tmp;
	}
    }
}
// hand-unroll innermost loop
void matmul_for_small_nk_2(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int kk,
			 int n)
{

    int i,j,k;
    double bcopy[n][kk];
    for(j=0;j<n;j++)
	    for(k=0;k<kk;k++)
		bcopy[j][k] = b[k][j];
    for(i=0;i<m;i++){
	double *ap=a[i];
	for(j=0;j<n;j++){
	    double *bp = bcopy[j];
	    double tmp=0.0;
	    for(k=0;k<kk;k+=8)
		tmp += ap[k]*bp[k]
		    + ap[k+1]*bp[k+1]
		    + ap[k+2]*bp[k+2]
		    + ap[k+3]*bp[k+3]
		    + ap[k+4]*bp[k+4]
		    + ap[k+5]*bp[k+5]
		    + ap[k+6]*bp[k+6]
		    + ap[k+7]*bp[k+7];
	    c[i][j]-=tmp;
	}
    }
}
// hand-unroll mid-loop
void matmul_for_small_nk_3(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int kk,
			 int n)
{

    int i,j,k;
    double bcopy[n][kk];
    for(j=0;j<n;j++)
	    for(k=0;k<kk;k++)
		bcopy[j][k] = b[k][j];
    for(i=0;i<m;i++){
	double *ap=a[i];
	for(j=0;j<n;j+=4){
	    double *bp = bcopy[j];
	    double *bpp = bcopy[j+1];
	    double *bp2 = bcopy[j+2];
	    double *bp3 = bcopy[j+3];
	    double tmp=0.0;
	    double tmp1=0.0;
	    double tmp2=0.0;
	    double tmp3=0.0;
	    for(k=0;k<kk;k+=8){
		tmp += ap[k]*bp[k]
		    + ap[k+1]*bp[k+1]
		    + ap[k+2]*bp[k+2]
		    + ap[k+3]*bp[k+3]
		    + ap[k+4]*bp[k+4]
		    + ap[k+5]*bp[k+5]
		    + ap[k+6]*bp[k+6]
		    + ap[k+7]*bp[k+7];
		tmp1 += ap[k]*bpp[k]
		    + ap[k+1]*bpp[k+1]
		    + ap[k+2]*bpp[k+2]
		    + ap[k+3]*bpp[k+3]
		    + ap[k+4]*bpp[k+4]
		    + ap[k+5]*bpp[k+5]
		    + ap[k+6]*bpp[k+6]
		    + ap[k+7]*bpp[k+7];
		tmp2 += ap[k]*bp2[k]
		    + ap[k+1]*bp2[k+1]
		    + ap[k+2]*bp2[k+2]
		    + ap[k+3]*bp2[k+3]
		    + ap[k+4]*bp2[k+4]
		    + ap[k+5]*bp2[k+5]
		    + ap[k+6]*bp2[k+6]
		    + ap[k+7]*bp2[k+7];
		tmp3 += ap[k]*bp3[k]
		    + ap[k+1]*bp3[k+1]
		    + ap[k+2]*bp3[k+2]
		    + ap[k+3]*bp3[k+3]
		    + ap[k+4]*bp3[k+4]
		    + ap[k+5]*bp3[k+5]
		    + ap[k+6]*bp3[k+6]
		    + ap[k+7]*bp3[k+7];
	    }
	    c[i][j]-=tmp;
	    c[i][j+1]-=tmp1;
	    c[i][j+2]-=tmp2;
	    c[i][j+3]-=tmp3;
	}
    }
}
// hand-unroll mid-loop by 2
void matmul_for_small_nk_4(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int kk,
			 int n)
{

    int i,j,k;
    double bcopy[n][kk];
    for(j=0;j<n;j++)
	    for(k=0;k<kk;k++)
		bcopy[j][k] = b[k][j];
    for(i=0;i<m;i++){
	double *ap=a[i];
	for(j=0;j<n;j+=2){
	    double *bp = bcopy[j];
	    double *bpp = bcopy[j+1];
	    double tmp=0.0;
	    double tmp1=0.0;
	    for(k=0;k<kk;k+=8){
		tmp += ap[k]*bp[k]
		    + ap[k+1]*bp[k+1]
		    + ap[k+2]*bp[k+2]
		    + ap[k+3]*bp[k+3]
		    + ap[k+4]*bp[k+4]
		    + ap[k+5]*bp[k+5]
		    + ap[k+6]*bp[k+6]
		    + ap[k+7]*bp[k+7];
		tmp1 += ap[k]*bpp[k]
		    + ap[k+1]*bpp[k+1]
		    + ap[k+2]*bpp[k+2]
		    + ap[k+3]*bpp[k+3]
		    + ap[k+4]*bpp[k+4]
		    + ap[k+5]*bpp[k+5]
		    + ap[k+6]*bpp[k+6]
		    + ap[k+7]*bpp[k+7];
	    }
	    c[i][j]-=tmp;
	    c[i][j+1]-=tmp1;
	}
    }
}
// use sse2 for dot product
void matmul_for_small_nk_5(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int kk,
			 int n)
{

    int i,j,k;
    int nh = n/2;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
    v2df acopy[kk];
    double acopyd[kk];
    for(j=0;j<nh;j++)
	    for(k=0;k<kk;k++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	double *ap=a[i];
	double *acp = (double*) acopy;
	register v2df tmp= (v2df){0.0,0.0};
	v2df * cp = (v2df*) (&(c[i][0]));
	for(k=0;k<kk;k+=4){
	     __builtin_prefetch((double*)a[i+4]+k,0);
	}
	for(j=0;j<n;j+=4){
	     __builtin_prefetch(c[i+4]+j,0);
	}

	for(k=0;k<kk;k+=2){
	    //	    v2df aa = *((v2df*)(ap+k));
	    //	    acopy[k]=__builtin_ia32_shufpd(aa,aa,0x0);
	    //	    acopy[k+1]= __builtin_ia32_shufpd(aa,aa,0x5);
	    acp[k*2]=acp[k*2+1]=ap[k];
	    acp[k*2+2]=acp[k*2+3]=ap[k+1];
	}
	for(j=0;j<nh;j++){
	    tmp = (v2df){0.0,0.0};
	    v2df * bp = bcopy2[j];
	    for(k=0;k<kk;k+=2){
		tmp += acopy[k]*bp[k]
		    +acopy[k+1]*bp[k+1]
#if 0
		    +acopy[k+2]*bp[k+2]
		    +acopy[k+3]*bp[k+3]
		    +acopy[k+4]*bp[k+4]
		    +acopy[k+5]*bp[k+5]
		    +acopy[k+6]*bp[k+6]
		    +acopy[k+7]*bp[k+7]
#endif
		    ;

	    }
	    cp[j] -= tmp;
	}
    }
}
// use sse2 for dot product
void matmul_for_small_nk_6(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int kk,
			 int n)
{

    int i,j,k;
    int nh = n/2;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
    v2df acopy[kk];
    v2df acopy2[kk];
    v2df acopy3[kk];
    v2df acopy4[kk];
    double *acp = (double*) acopy;
    double *acp2 = (double*) acopy2;
    double *acp3 = (double*) acopy3;
    double *acp4 = (double*) acopy4;
    for(j=0;j<nh;j++)
	    for(k=0;k<kk;k++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i+=4){
	double *ap=a[i];
	double *ap2=a[i+1];
	double *ap3=a[i+2];
	double *ap4=a[i+3];
	register v2df tmp, tmp2, tmp3, tmp4;
	v2df * cp = (v2df*) (&(c[i][0]));
	v2df * cp2 = (v2df*) (&(c[i+1][0]));
	v2df * cp3 = (v2df*) (&(c[i+2][0]));
	v2df * cp4 = (v2df*) (&(c[i+3][0]));
	for(k=0;k<kk;k+=4){
	     __builtin_prefetch((double*)a[i+4]+k,0);
	     __builtin_prefetch((double*)a[i+5]+k,0);
	     __builtin_prefetch((double*)a[i+6]+k,0);
	     __builtin_prefetch((double*)a[i+7]+k,0);
	     __builtin_prefetch(c[i+4]+j,0);
	     __builtin_prefetch(c[i+5]+j,0);
	     __builtin_prefetch(c[i+6]+j,0);
	     __builtin_prefetch(c[i+7]+j,0);
	}

	for(k=0;k<kk;k+=2){
	    v2df * aa = (v2df*)(ap+k);
	    acopy [k]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	    acopy[k+1]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	    aa = (v2df*)(ap2+k);
	    acopy2 [k]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	    acopy2[k+1]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	    //   acp[k*2]=acp[k*2+1]=ap[k];
	    //	    acp[k*2+2]=acp[k*2+3]=ap[k+1];
	    //	    acp2[k*2]=acp2[k*2+1]=ap2[k];
	    //	    acp2[k*2+2]=acp2[k*2+3]=ap2[k+1];
	}
	for(k=0;k<kk;k+=2){
	    acp3[k*2]=acp3[k*2+1]=ap3[k];
	    acp3[k*2+2]=acp3[k*2+3]=ap3[k+1];
	    acp4[k*2]=acp4[k*2+1]=ap4[k];
	    acp4[k*2+2]=acp4[k*2+3]=ap4[k+1];
	}
	for(j=0;j<nh;j++){
	    tmp = tmp2= tmp3= tmp4= (v2df){0.0,0.0};
	    v2df * bp = bcopy2[j];
#if 0
	    for(k=0;k<kk;k+=4){
		tmp += acopy[k]*bp[k]
		    +acopy[k+1]*bp[k+1]
		    +acopy[k+2]*bp[k+2]
		    +acopy[k+3]*bp[k+3];
		tmp2 += acopy2[k]*bp[k]
		    +acopy2[k+1]*bp[k+1]
		    +acopy2[k+2]*bp[k+2]
		    +acopy2[k+3]*bp[k+3];
		tmp3 += acopy3[k]*bp[k]
		    +acopy3[k+1]*bp[k+1]
		    +acopy3[k+2]*bp[k+2]
		    +acopy3[k+3]*bp[k+3];
		tmp4 += acopy4[k]*bp[k]
		    +acopy4[k+1]*bp[k+1]
		    +acopy4[k+2]*bp[k+2]
		    +acopy4[k+3]*bp[k+3];

	    }
#endif
	    for(k=0;k<kk;k+=2){
		tmp += acopy[k]*bp[k]
		    +acopy[k+1]*bp[k+1];
		tmp2 += acopy2[k]*bp[k]
		    +acopy2[k+1]*bp[k+1];
		tmp3 += acopy3[k]*bp[k]
		    +acopy3[k+1]*bp[k+1];
		tmp4 += acopy4[k]*bp[k]
		    +acopy4[k+1]*bp[k+1];

	    }
	    cp[j] -= tmp;
	    cp2[j] -= tmp2;
	    cp3[j] -= tmp3;
	    cp4[j] -= tmp4;
	}
    }
}
void matmul_for_small_nk7(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int kk,
			 int n)
{

    int i,j;
    int nh = n/2;
    register int k;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
    v2df acopy[kk];
    v2df acopy2[kk];
    double *acp = (double*) acopy;
    double *acp2 = (double*) acopy2;
    unsigned long bpcount, apcount, dotcount;
    bpcount= apcount= dotcount=0;
    //    BEGIN_TSC;
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    //    END_TSC(bpcount);
    for(i=0;i<m;i+=2){
	//	BEGIN_TSC;
	double *ap=a[i];
	double *ap2=a[i+1];
	register v2df tmp, tmp2;
	v2df * cp = (v2df*) (&(c[i][0]));
	v2df * cp2 = (v2df*) (&(c[i+1][0]));

	for(k=0;k<kk;k+=2){
	    v2df * aa = (v2df*)(ap+k);
	     __builtin_prefetch((double*)a[i+4]+k,0);
	     __builtin_prefetch((double*)a[i+5]+k,0);
	    acopy [k]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	    acopy[k+1]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	    aa = (v2df*)(ap2+k);
	    acopy2 [k]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	    acopy2[k+1]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	    //   acp[k*2]=acp[k*2+1]=ap[k];
	    //	    acp[k*2+2]=acp[k*2+3]=ap[k+1];
	    //	    acp2[k*2]=acp2[k*2+1]=ap2[k];
	    //	    acp2[k*2+2]=acp2[k*2+3]=ap2[k+1];
	}
	//	END_TSC(apcount);
	//	BEGIN_TSC;
	for(j=0;j<nh;j++){
	    tmp = tmp2=  (v2df){0.0,0.0};
	    v2df ctmp= cp[j];
	    v2df ctmp2 = cp2[j] ;

	    v2df * bp = bcopy2[j];
	    __builtin_prefetch(c[i+4]+j,0);
	    __builtin_prefetch(c[i+5]+j,0);
	    for(k=0;k<kk;k+=8){
		int k2 = k+4;
		v2df *avp = acopy+k;
		v2df *avp2 = acopy2+k;
		v2df *bvp = bp+k;
		tmp += avp[0]*bvp[0];
		tmp2 += avp2[0]*bvp[0];
		tmp +=avp[1]*bvp[1];
		tmp2+=avp2[1]*bvp[1];
		tmp +=avp[2]*bvp[2];
		tmp2+=avp2[2]*bvp[2];
		tmp +=avp[3]*bvp[3];
		tmp2+=avp2[3]*bvp[3];
		tmp += avp[4]*bvp[4];
		tmp2 += avp2[4]*bvp[4];
		tmp +=avp[5]*bvp[5];
		tmp2+=avp2[5]*bvp[5];
		tmp +=avp[6]*bvp[6];
		tmp2+=avp2[6]*bvp[6];
		tmp +=avp[7]*bvp[7];
		tmp2+=avp2[7]*bvp[7];
	    }

#if 0
	    for(k=0;k<kk;k+=8){
		int k2 = k+4;
		tmp += acopy[k]*bp[k];
		tmp2 += acopy2[k]*bp[k];
		tmp +=acopy[k+1]*bp[k+1];
		tmp2+=acopy2[k+1]*bp[k+1];
		tmp +=acopy[k+2]*bp[k+2];
		tmp2+=acopy2[k+2]*bp[k+2];
		tmp +=acopy[k+3]*bp[k+3];
		tmp2+=acopy2[k+3]*bp[k+3];
		tmp += acopy[k2]*bp[k2];
		tmp2 += acopy2[k2]*bp[k2];
		tmp +=acopy[k2+1]*bp[k2+1];
		tmp2+=acopy2[k2+1]*bp[k2+1];
		tmp +=acopy[k2+2]*bp[k2+2];
		tmp2+=acopy2[k2+2]*bp[k2+2];
		tmp +=acopy[k2+3]*bp[k2+3];
		tmp2+=acopy2[k2+3]*bp[k2+3];
	    }
#endif
#if 0
	    for(k=0;k<kk;k+=2){
		tmp += acopy[k]*bp[k];
		tmp 	    __builtin_prefetch(c[i+4+(j&1)]+j,0);
+=acopy[k+1]*bp[k+1];
		tmp2 += acopy2[k]*bp[k];
		tmp2+=acopy2[k+1]*bp[k+1];
	    }
#endif
	    cp[j] = ctmp -tmp;
	    cp2[j] = ctmp2 -tmp2;

	}
	//	END_TSC(dotcount);

    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //	   (double)bpcount, (double)apcount, (double)dotcount);

}

// XMM registers
#define X0	"%xmm0"
#define X1	"%xmm1"
#define X2	"%xmm2"
#define X3	"%xmm3"
#define X4	"%xmm4"
#define X5	"%xmm5"
#define X6	"%xmm6"
#define X7	"%xmm7"
#define X8		"%xmm8"
#define X9		"%xmm9"
#define X10		"%xmm10"
#define X11		"%xmm11"
#define X12		"%xmm12"
#define X13		"%xmm13"
#define X14		"%xmm14"
#define X15		"%xmm15"

#define LOADPD(mem, reg) asm("movapd %0, %"reg::"m"(mem));
#define STORPD(reg, mem) asm("movapd %"reg " , %0"::"m"(mem));
#define MOVNTPD(reg, mem) asm("movntpd %"reg " , %0"::"m"(mem));
#define MOVAPD(src, dst) asm("movapd " src "," dst);
#define MOVQ(src, dst) asm("movq " src "," dst);
#define BCAST0(reg) asm("shufpd $0x00, " reg ","  reg);
#define BCAST1(reg) asm("shufpd $0xff, " reg ","  reg);
#define MULPD(src, dst) asm("mulpd " src "," dst);
#define ADDPD(src, dst) asm("addpd " src ","  dst);
#define SUBPD(src, dst) asm("subpd "  src "," dst);

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

    int i,j;
    int kk = 8;
    int nh = n/2;
    register int k;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
#define PREFETCHL 32
    for(i=0;i<PREFETCHL;i++){
	__builtin_prefetch((double*)a[i],0,0);
	__builtin_prefetch(c[i+8],1,0);
    }

    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	//	BEGIN_TSC;
	//	v2df acopy[8];
	v2df *ap = (v2df*) a[i];
	v2df * cp = (v2df*) (&(c[i][0]));

	__builtin_prefetch((double*)a[i+PREFETCHL],0,0);

	int k;

	for(j=0;j<nh;j+=4){
	    __builtin_prefetch(c[i+PREFETCHL]+j,1,3);
	    v2df * bvp0 = bcopy2[j];
	    v2df * bvp1 = bcopy2[j+1];
	    v2df * bvp2 = bcopy2[j+2];
	    v2df * bvp3 = bcopy2[j+3];
	    LOADPD(cp[j],X12);
	    LOADPD(cp[j+1],X13);
	    LOADPD(cp[j+2],X14);
	    LOADPD(cp[j+3],X15);

	    LOADPD(ap[0],X0);
	    LOADPD(bvp0[0],X4);
	    LOADPD(bvp1[0],X5);
	    LOADPD(bvp2[0],X6);
	    LOADPD(bvp3[0],X7);
	    LOADPD(bvp0[1],X8);
	    LOADPD(bvp1[1],X9);
	    LOADPD(bvp2[1],X10);
	    LOADPD(bvp3[1],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);
	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[1],X0);

	    LOADPD(bvp0[2],X4);
	    LOADPD(bvp1[2],X5);
	    LOADPD(bvp2[2],X6);
	    LOADPD(bvp3[2],X7);
	    LOADPD(bvp0[3],X8);
	    LOADPD(bvp1[3],X9);
	    LOADPD(bvp2[3],X10);
	    LOADPD(bvp3[3],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);
	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[2],X0);
	    LOADPD(bvp0[4],X4);
	    LOADPD(bvp1[4],X5);
	    LOADPD(bvp2[4],X6);
	    LOADPD(bvp3[4],X7);
	    LOADPD(bvp0[5],X8);
	    LOADPD(bvp1[5],X9);
	    LOADPD(bvp2[5],X10);
	    LOADPD(bvp3[5],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);
	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[3],X0);

	    LOADPD(bvp0[6],X4);
	    LOADPD(bvp1[6],X5);
	    LOADPD(bvp2[6],X6);
	    LOADPD(bvp3[6],X7);
	    LOADPD(bvp0[7],X8);
	    LOADPD(bvp1[7],X9);
	    LOADPD(bvp2[7],X10);
	    LOADPD(bvp3[7],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);
	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);
	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);
	    STORPD(X12,cp[j+0]);
	    STORPD(X13,cp[j+1]);
	    STORPD(X14,cp[j+2]);
	    STORPD(X15,cp[j+3]);
	}
    }
}

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

    int i,j;
    int kk = 16;
    int nh = n/2;
    register int k;
    v2df bcopy2[nh][kk];
#ifdef     PREFETCHL
#undef PREFETCHL
#endif
#define PREFETCHL 32

    for(i=0;i<PREFETCHL;i++){
	__builtin_prefetch((double*)a[i],0,0);
	__builtin_prefetch((double*)a[i]+8,0,0);
	__builtin_prefetch(c[i+8],1,0);
	__builtin_prefetch(c[i+8]+8,1,0);
    }
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i+=2){
	//	BEGIN_TSC;
	v2df *ap = (v2df*) a[i];
	v2df * cp = (v2df*) (&(c[i][0]));
	v2df *app = (v2df*) a[i+1];
	v2df * cpp = (v2df*) (&(c[i+1][0]));

	__builtin_prefetch((double*)a[i+PREFETCHL],0,0);
	__builtin_prefetch((double*)a[i+PREFETCHL]+8,0,0);
	__builtin_prefetch((double*)a[i+PREFETCHL+1],0,0);
	__builtin_prefetch((double*)a[i+PREFETCHL+1]+8,0,0);

	int k;

	for(j=0;j<nh;j+=2){
	    __builtin_prefetch(c[i+PREFETCHL]+j,1,0);
	    v2df * bvp0 = bcopy2[j];
	    v2df * bvp1 = bcopy2[j+1];
	    LOADPD(cp[j],X12);
	    LOADPD(cp[j+1],X13);
	    LOADPD(cpp[j],X14);
	    LOADPD(cpp[j+1],X15);

	    for(k=0;k<8;k++){

		LOADPD(ap[k],X0);
		LOADPD(app[k],X2);
		MOVAPD(X0,X1);
		BCAST0(X0);
		BCAST1(X1);
		MOVAPD(X2,X3);
		BCAST0(X2);
		BCAST1(X3);
		LOADPD(bvp0[k*2],X4);
		MOVAPD(X4,X6);
		MULPD(X0,X4);
		SUBPD(X4,X12);
		LOADPD(bvp1[k*2],X5);
		MOVAPD(X5,X7);
		MULPD(X0,X5);
		SUBPD(X5,X13);
		LOADPD(bvp0[k*2+1],X8);
		MOVAPD(X8,X10);
		MULPD(X1,X8);
		SUBPD(X8,X12);
		LOADPD(bvp1[k*2+1],X9);
		MOVAPD(X9,X11);
		MULPD(X1,X9);
		SUBPD(X9,X13);
		MULPD(X2,X6);
		SUBPD(X6,X14);
		MULPD(X2,X7);
		SUBPD(X7,X15);
		MULPD(X3,X10);
		SUBPD(X10,X14);
		MULPD(X3,X11);
		SUBPD(X11,X15);
	    }

	    STORPD(X12,cp[j+0]);
	    STORPD(X13,cp[j+1]);
	    STORPD(X14,cpp[j+0]);
	    STORPD(X15,cpp[j+1]);
	}
    }
}
void matmul_for_nk16_0c(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j;
    int kk = 16;
    int nh = n/2;
    register int k;
    v2df bcopy2[nh][kk];
#ifdef     PREFETCHL
#undef PREFETCHL
#endif
#define PREFETCHL 16

    for(i=0;i<PREFETCHL;i++){
	__builtin_prefetch((double*)a[i],0,0);
	__builtin_prefetch((double*)a[i]+8,0,0);
	__builtin_prefetch(c[i+8],1,0);
	__builtin_prefetch(c[i+8]+8,1,0);
    }
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	//	BEGIN_TSC;
	v2df *ap = (v2df*) a[i];
	v2df * cp = (v2df*) (&(c[i][0]));

	__builtin_prefetch((double*)a[i+PREFETCHL],0,0);
	__builtin_prefetch((double*)a[i+PREFETCHL]+8,0,0);

	int k;

	for(j=0;j<nh;j+=4){
	    __builtin_prefetch(c[i+PREFETCHL]+j,1,0);
	    v2df * bvp0 = bcopy2[j];
	    v2df * bvp1 = bcopy2[j+1];
	    v2df * bvp2 = bcopy2[j+2];
	    v2df * bvp3 = bcopy2[j+3];
	    LOADPD(cp[j],X12);
	    LOADPD(cp[j+1],X13);
	    LOADPD(cp[j+2],X14);
	    LOADPD(cp[j+3],X15);

	    for(k=0;k<8;k++){

		LOADPD(ap[k],X0);
		LOADPD(bvp0[k*2],X4);
		LOADPD(bvp1[k*2],X5);
		LOADPD(bvp2[k*2],X6);
		LOADPD(bvp3[k*2],X7);
		MOVAPD(X0,X1);
		BCAST0(X0);
		BCAST1(X1);
		LOADPD(bvp0[k*2+1],X8);
		LOADPD(bvp1[k*2+1],X9);
		LOADPD(bvp2[k*2+1],X10);
		LOADPD(bvp3[k*2+1],X11);
		MULPD(X0,X4);
		MULPD(X0,X5);
		MULPD(X0,X6);
		MULPD(X0,X7);
		MULPD(X1,X8);
		MULPD(X1,X9);
		MULPD(X1,X10);
		MULPD(X1,X11);
		SUBPD(X4,X12);
		SUBPD(X5,X13);
		SUBPD(X6,X14);
		SUBPD(X7,X15);
		SUBPD(X8,X12);
		SUBPD(X9,X13);
		SUBPD(X10,X14);
		SUBPD(X11,X15);
	    }

	    STORPD(X12,cp[j+0]);
	    STORPD(X13,cp[j+1]);
	    STORPD(X14,cp[j+2]);
	    STORPD(X15,cp[j+3]);
	}
    }
}
void matmul_for_nk32_0(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j;
    int kk = 32;
    int nh = n/2;
    register int k;
    v2df bcopy2[nh][kk];
#ifdef     PREFETCHL
#undef PREFETCHL
#endif
#define PREFETCHL 8

    for(i=0;i<PREFETCHL;i++){
	__builtin_prefetch((double*)a[i],0,0);
	__builtin_prefetch((double*)a[i]+8,0,0);
	__builtin_prefetch((double*)a[i]+16,0,0);
	__builtin_prefetch((double*)a[i]+24,0,0);
	__builtin_prefetch(c[i+8],1,0);
	__builtin_prefetch(c[i+8]+8,1,0);
	__builtin_prefetch(c[i+8]+16,1,0);
	__builtin_prefetch(c[i+8]+24,1,0);
    }
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	//	BEGIN_TSC;
	v2df *ap = (v2df*) a[i];
	v2df * cp = (v2df*) (&(c[i][0]));

	__builtin_prefetch((double*)a[i+PREFETCHL],0,0);
	__builtin_prefetch((double*)a[i+PREFETCHL]+8,0,0);
	__builtin_prefetch((double*)a[i+PREFETCHL]+16,0,0);
	__builtin_prefetch((double*)a[i+PREFETCHL]+24,0,0);

	int k;

	for(j=0;j<nh;j+=4){
	    __builtin_prefetch(c[i+PREFETCHL]+j,1,0);
	    v2df * bvp0 = bcopy2[j];
	    v2df * bvp1 = bcopy2[j+1];
	    v2df * bvp2 = bcopy2[j+2];
	    v2df * bvp3 = bcopy2[j+3];
	    LOADPD(cp[j],X12);
	    LOADPD(cp[j+1],X13);
	    LOADPD(cp[j+2],X14);
	    LOADPD(cp[j+3],X15);

	    for(k=0;k<16;k++){

		LOADPD(ap[k],X0);
		LOADPD(bvp0[k*2],X4);
		LOADPD(bvp1[k*2],X5);
		LOADPD(bvp2[k*2],X6);
		LOADPD(bvp3[k*2],X7);
		MOVAPD(X0,X1);
		BCAST0(X0);
		BCAST1(X1);
		MULPD(X0,X4);
		MULPD(X0,X5);
		MULPD(X0,X6);
		MULPD(X0,X7);
		LOADPD(bvp0[k*2+1],X8);
		LOADPD(bvp1[k*2+1],X9);
		LOADPD(bvp2[k*2+1],X10);
		LOADPD(bvp3[k*2+1],X11);
		MULPD(X1,X8);
		MULPD(X1,X9);
		MULPD(X1,X10);
		MULPD(X1,X11);
		SUBPD(X4,X12);
		SUBPD(X5,X13);
		SUBPD(X6,X14);
		SUBPD(X7,X15);
		SUBPD(X8,X12);
		SUBPD(X9,X13);
		SUBPD(X10,X14);
		SUBPD(X11,X15);
	    }

	    STORPD(X12,cp[j+0]);
	    STORPD(X13,cp[j+1]);
	    STORPD(X14,cp[j+2]);
	    STORPD(X15,cp[j+3]);
	}
    }
}

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

    int i,j;
    int kk = 16;
    int nh = n/2;
    register int k;
    v2df bcopy2[nh][kk];
#ifdef     PREFETCHL
#undef PREFETCHL
#endif
#define PREFETCHL 16

    for(i=0;i<PREFETCHL;i++){
	__builtin_prefetch((double*)a[i],0,0);
	__builtin_prefetch((double*)a[i]+8,0,0);
	__builtin_prefetch(c[i+8],1,0);
	__builtin_prefetch(c[i+8]+8,1,0);
    }
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	//	BEGIN_TSC;
	v2df *ap = (v2df*) a[i];
	v2df * cp = (v2df*) (&(c[i][0]));

	__builtin_prefetch((double*)a[i+PREFETCHL],0,0);
	__builtin_prefetch((double*)a[i+PREFETCHL]+8,0,0);

	int k;

	for(j=0;j<nh;j+=4){
	    __builtin_prefetch(c[i+PREFETCHL]+j,1,0);
	    v2df * bvp0 = bcopy2[j];
	    v2df * bvp1 = bcopy2[j+1];
	    v2df * bvp2 = bcopy2[j+2];
	    v2df * bvp3 = bcopy2[j+3];
	    LOADPD(cp[j],X12);
	    LOADPD(cp[j+1],X13);
	    LOADPD(cp[j+2],X14);
	    LOADPD(cp[j+3],X15);

	    LOADPD(ap[0],X0);
	    LOADPD(bvp0[0],X4);
	    LOADPD(bvp1[0],X5);
	    LOADPD(bvp2[0],X6);
	    LOADPD(bvp3[0],X7);
	    LOADPD(bvp0[1],X8);
	    LOADPD(bvp1[1],X9);
	    LOADPD(bvp2[1],X10);
	    LOADPD(bvp3[1],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[1],X0);

	    LOADPD(bvp0[2],X4);
	    LOADPD(bvp1[2],X5);
	    LOADPD(bvp2[2],X6);
	    LOADPD(bvp3[2],X7);
	    LOADPD(bvp0[3],X8);
	    LOADPD(bvp1[3],X9);
	    LOADPD(bvp2[3],X10);
	    LOADPD(bvp3[3],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[2],X0);
	    LOADPD(bvp0[4],X4);
	    LOADPD(bvp1[4],X5);
	    LOADPD(bvp2[4],X6);
	    LOADPD(bvp3[4],X7);
	    LOADPD(bvp0[5],X8);
	    LOADPD(bvp1[5],X9);
	    LOADPD(bvp2[5],X10);
	    LOADPD(bvp3[5],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[3],X0);
	    LOADPD(bvp0[6],X4);
	    LOADPD(bvp1[6],X5);
	    LOADPD(bvp2[6],X6);
	    LOADPD(bvp3[6],X7);
	    LOADPD(bvp0[7],X8);
	    LOADPD(bvp1[7],X9);
	    LOADPD(bvp2[7],X10);
	    LOADPD(bvp3[7],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[4],X0);

	    LOADPD(bvp0[8],X4);
	    LOADPD(bvp1[8],X5);
	    LOADPD(bvp2[8],X6);
	    LOADPD(bvp3[8],X7);
	    LOADPD(bvp0[9],X8);
	    LOADPD(bvp1[9],X9);
	    LOADPD(bvp2[9],X10);
	    LOADPD(bvp3[9],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);
	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[5],X0);
	    LOADPD(bvp0[10],X4);
	    LOADPD(bvp1[10],X5);
	    LOADPD(bvp2[10],X6);
	    LOADPD(bvp3[10],X7);
	    LOADPD(bvp0[11],X8);
	    LOADPD(bvp1[11],X9);
	    LOADPD(bvp2[11],X10);
	    LOADPD(bvp3[11],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[6],X0);
	    LOADPD(bvp0[12],X4);
	    LOADPD(bvp1[12],X5);
	    LOADPD(bvp2[12],X6);
	    LOADPD(bvp3[12],X7);
	    LOADPD(bvp0[13],X8);
	    LOADPD(bvp1[13],X9);
	    LOADPD(bvp2[13],X10);
	    LOADPD(bvp3[13],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    LOADPD(ap[7],X0);
	    LOADPD(bvp0[14],X4);
	    LOADPD(bvp1[14],X5);
	    LOADPD(bvp2[14],X6);
	    LOADPD(bvp3[14],X7);
	    LOADPD(bvp0[15],X8);
	    LOADPD(bvp1[15],X9);
	    LOADPD(bvp2[15],X10);
	    LOADPD(bvp3[15],X11);
	    MOVAPD(X0,X1);
	    BCAST0(X0);
	    BCAST1(X1);

	    MULPD(X0,X4);
	    MULPD(X0,X5);
	    MULPD(X0,X6);
	    MULPD(X0,X7);
	    MULPD(X1,X8);
	    MULPD(X1,X9);
	    MULPD(X1,X10);
	    MULPD(X1,X11);
	    SUBPD(X4,X12);
	    SUBPD(X5,X13);
	    SUBPD(X6,X14);
	    SUBPD(X7,X15);

	    SUBPD(X8,X12);
	    SUBPD(X9,X13);
	    SUBPD(X10,X14);
	    SUBPD(X11,X15);

	    STORPD(X12,cp[j+0]);
	    STORPD(X13,cp[j+1]);
	    STORPD(X14,cp[j+2]);
	    STORPD(X15,cp[j+3]);
	}
    }
}

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

    int i,j;
    int kk = 8;
    int nh = n/2;
    register int k;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
    v2df acopy[kk];
    v2df acopy2[kk];
    v2df awork[4];
    v2df awork2[4];
    double *acp = (double*) acopy;
    double *acp2 = (double*) acopy2;
    unsigned long bpcount, apcount, dotcount;
    bpcount= apcount= dotcount=0;
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	//	BEGIN_TSC;
	double *ap=a[i];
	register v2df tmp, tmp2;
	v2df * cp = (v2df*) (&(c[i][0]));

	v2df * aa = (v2df*)(ap);
	__builtin_prefetch((double*)a[i+8],0,0);

	v2df acopy0=(v2df){a[i][0], a[i][0]};
	v2df acopy1=(v2df){a[i][1], a[i][1]};
	v2df acopy2=(v2df){a[i][2], a[i][2]};
	v2df acopy3=(v2df){a[i][3], a[i][3]};
	v2df acopy4=(v2df){a[i][4], a[i][4]};
	v2df acopy5=(v2df){a[i][5], a[i][5]};
	v2df acopy6=(v2df){a[i][6], a[i][6]};
	v2df acopy7=(v2df){a[i][7], a[i][7]};
	v2df zero=(v2df){0.0, 0.0};

	LOADPD(acopy0,X0);
	LOADPD(acopy1,X1);
	LOADPD(acopy2,X2);
	LOADPD(acopy3,X3);
	LOADPD(acopy4,X4);
	LOADPD(acopy5,X5);
	LOADPD(acopy6,X6);
	LOADPD(acopy7,X7);

	for(j=0;j<nh;j++){
	    __builtin_prefetch(c[i+8]+j,1,0);
	    v2df * bvp = bcopy2[j];
	    LOADPD(cp[j],X14);
	    LOADPD(bvp[0],X8);
	    LOADPD(bvp[1],X9);
	    MULPD(X0,X8);
	    MULPD(X1,X9);
	    LOADPD(bvp[2],X10);
	    LOADPD(bvp[3],X11);
	    ADDPD(X9,X8);
	    MULPD(X2,X10);
	    MULPD(X3,X11);
	    ADDPD(X11,X10);
	    LOADPD(bvp[4],X9);
	    LOADPD(bvp[5],X11);
	    LOADPD(bvp[6],X12);
	    LOADPD(bvp[7],X13);
	    MULPD(X4,X9);
	    MULPD(X5,X11);
	    ADDPD(X10,X8);
	    ADDPD(X11,X9);
	    MULPD(X6,X12);
	    MULPD(X7,X13);
	    ADDPD(X13,X12);
	    ADDPD(X9,X8);
	    ADDPD(X12,X8);
	    SUBPD(X8,X14);
	    STORPD(X14,cp[j]);
	}
    }
}

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

    int i,j;
    int kk = 8;
    int nh = n/2;
    register int k;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
    v2df acopy[kk];
    v2df acopy2[kk];
    v2df awork[4];
    v2df awork2[4];
    double *acp = (double*) acopy;
    double *acp2 = (double*) acopy2;
    unsigned long bpcount, apcount, dotcount;
    bpcount= apcount= dotcount=0;
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	//	BEGIN_TSC;
	double *ap=a[i];
	register v2df tmp, tmp2;
	v2df * cp = (v2df*) (&(c[i][0]));

	v2df * aa = (v2df*)(ap);
	__builtin_prefetch((double*)a[i+8],0,0);

	register v2df acopy0=(v2df){a[i][0], a[i][0]};
	register v2df acopy1=(v2df){a[i][1], a[i][1]};
	register v2df acopy2=(v2df){a[i][2], a[i][2]};
	register v2df acopy3=(v2df){a[i][3], a[i][3]};
	register v2df acopy4=(v2df){a[i][4], a[i][4]};
	register v2df acopy5=(v2df){a[i][5], a[i][5]};
	register v2df acopy6=(v2df){a[i][6], a[i][6]};
	register v2df acopy7=(v2df){a[i][7], a[i][7]};

	for(j=0;j<nh;j++){
	    tmp =  (v2df){0.0,0.0};
	    v2df ctmp= cp[j];

	    v2df * bp = bcopy2[j];
	    __builtin_prefetch(c[i+4]+j,1,0);
	    v2df *bvp = bp;
	    tmp += acopy0*bvp[0];
	    tmp +=acopy1*bvp[1];
	    tmp +=acopy2*bvp[2];
	    tmp +=acopy3*bvp[3];
	    tmp +=acopy4*bvp[4];
	    tmp +=acopy5*bvp[5];
	    tmp +=acopy6*bvp[6];
	    tmp +=acopy7*bvp[7];
	    cp[j] = ctmp -tmp;
	}
    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //	   (double)bpcount, (double)apcount, (double)dotcount);

}
void matmul_for_nk8_0b(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j;
    int kk = 8;
    int nh = n/2;
    register int k;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
    v2df acopy[kk];
    v2df acopy2[kk];
    v2df awork[4];
    v2df awork2[4];
    double *acp = (double*) acopy;
    double *acp2 = (double*) acopy2;
    unsigned long bpcount, apcount, dotcount;
    bpcount= apcount= dotcount=0;
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i++){
	//	BEGIN_TSC;
	double *ap=a[i];
	register v2df tmp, tmp2;
	v2df * cp = (v2df*) (&(c[i][0]));

	v2df * aa = (v2df*)(ap);
	__builtin_prefetch((double*)a[i+8],0,0);

	acopy[0]=(v2df){a[i][0], a[i][0]};
	acopy[1]=(v2df){a[i][1], a[i][1]};
	acopy[2]=(v2df){a[i][2], a[i][2]};
	acopy[3]=(v2df){a[i][3], a[i][3]};
	acopy[4]=(v2df){a[i][4], a[i][4]};
	acopy[5]=(v2df){a[i][5], a[i][5]};
	acopy[6]=(v2df){a[i][6], a[i][6]};
	acopy[7]=(v2df){a[i][7], a[i][7]};

	for(j=0;j<nh;j++){
	    tmp = tmp2=  (v2df){0.0,0.0};
	    v2df ctmp= cp[j];

	    v2df * bp = bcopy2[j];
	    __builtin_prefetch(c[i+4]+j,1,0);
	    v2df *avp = acopy;
	    v2df *bvp = bp;
	    tmp += avp[0]*bvp[0];
	    tmp +=avp[1]*bvp[1];
	    tmp +=avp[2]*bvp[2];
	    tmp +=avp[3]*bvp[3];
	    tmp += avp[4]*bvp[4];
	    tmp +=avp[5]*bvp[5];
	    tmp +=avp[6]*bvp[6];
	    tmp +=avp[7]*bvp[7];
	    cp[j] = ctmp -tmp;
	}
    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //	   (double)bpcount, (double)apcount, (double)dotcount);

}
void matmul_for_nk8_0a(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j;
    int kk = 8;
    int nh = n/2;
    register int k;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
    v2df acopy[kk];
    v2df acopy2[kk];
    v2df awork[4];
    v2df awork2[4];
    double *acp = (double*) acopy;
    double *acp2 = (double*) acopy2;
    unsigned long bpcount, apcount, dotcount;
    bpcount= apcount= dotcount=0;
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    for(i=0;i<m;i+=2){
	//	BEGIN_TSC;
	double *ap=a[i];
	double *ap2=a[i+1];
	register v2df tmp, tmp2;
	v2df * cp = (v2df*) (&(c[i][0]));
	v2df * cp2 = (v2df*) (&(c[i+1][0]));

	v2df * aa = (v2df*)(ap);
	__builtin_prefetch((double*)a[i+8],0,0);
	__builtin_prefetch((double*)a[i+9],0,0);

	acopy[0]=(v2df){a[i][0], a[i][0]};
	acopy[1]=(v2df){a[i][1], a[i][1]};
	acopy[2]=(v2df){a[i][2], a[i][2]};
	acopy[3]=(v2df){a[i][3], a[i][3]};
	acopy[4]=(v2df){a[i][4], a[i][4]};
	acopy[5]=(v2df){a[i][5], a[i][5]};
	acopy[6]=(v2df){a[i][6], a[i][6]};
	acopy[7]=(v2df){a[i][7], a[i][7]};

	aa = (v2df*)(ap2);
	acopy2[0]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	acopy2[1]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	aa++;
	acopy2[2]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	acopy2[3]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	aa++;
	acopy2[4]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	acopy2[5]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	aa++;
	acopy2[6]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	acopy2[7]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	aa++;
	for(j=0;j<nh;j++){
	    tmp = tmp2=  (v2df){0.0,0.0};
	    v2df ctmp= cp[j];
	    v2df ctmp2 = cp2[j] ;

	    v2df * bp = bcopy2[j];
	    __builtin_prefetch(c[i+4]+j,1,0);
	    __builtin_prefetch(c[i+5]+j,1,0);
	    v2df *avp = acopy;
	    v2df *avp2 = acopy2;
	    v2df *bvp = bp;
	    tmp += avp[0]*bvp[0];
	    tmp2 += avp2[0]*bvp[0];
	    tmp +=avp[1]*bvp[1];
	    tmp2+=avp2[1]*bvp[1];
	    tmp +=avp[2]*bvp[2];
	    tmp2+=avp2[2]*bvp[2];
	    tmp +=avp[3]*bvp[3];
	    tmp2+=avp2[3]*bvp[3];
	    tmp += avp[4]*bvp[4];
	    tmp2 += avp2[4]*bvp[4];
	    tmp +=avp[5]*bvp[5];
	    tmp2+=avp2[5]*bvp[5];
	    tmp +=avp[6]*bvp[6];
	    tmp2+=avp2[6]*bvp[6];
	    tmp +=avp[7]*bvp[7];
	    tmp2+=avp2[7]*bvp[7];
	    cp[j] = ctmp -tmp;
	    cp2[j] = ctmp2 -tmp2;

	}

    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //	   (double)bpcount, (double)apcount, (double)dotcount);

}
void matmul_for_nk8_1(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j,k;
    const int kk = 8;
    const int kh = kk/2;
    int nh = n/2;
    v2df bcopy[n][kh];
    v2df acopy[kk][kh];
    unsigned long bpcount, apcount, dotcount;
    bpcount= apcount= dotcount=0;
    for(j=0;j<n;j++)
	for(k=0;k<kh;k++)
	    bcopy[j][k] = (v2df){b[k*2][j],b[k*2+1][j]};
    //    printf("copy b end\n");
    for(i=0;i<m;i+=kk){
	for(k=0;k<kk;k++){
	    v2df  *ak = (v2df*)(a[i+k]);
	    v2df * awp =acopy+k;
	    awp[0]=ak[0];
	    awp[1]=ak[1];
	    awp[2]=ak[2];
	    awp[3]=ak[3];
	}
	//	printf("copy a end\n");
	for(k=0;k<kk;k++){
	    v2u tmp, tmp1;
	    v2df * ap = acopy[k];
	    for(j=0;j<n;j+=2){
		tmp.v = ap[0]*bcopy[j][0]
		    + ap[1]*bcopy[j][1]
		    + ap[2]*bcopy[j][2]
		    + ap[3]*bcopy[j][3];
		tmp1.v = ap[0]*bcopy[j+1][0]
		    + ap[1]*bcopy[j+1][1]
		    + ap[2]*bcopy[j+1][2]
		    + ap[3]*bcopy[j+1][3];
		c[k+i][j] -= tmp.s[0]+tmp.s[1];
		c[k+i][j+1] -= tmp1.s[0]+tmp1.s[1];
	    }
	}
	//	printf("calc c end\n");
    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //	   (double)bpcount, (double)apcount, (double)dotcount);
}
void matmul_for_nk8_2(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j,k;
    const int kk = 8;
    const int kh = kk/2;
    int nh = n/2;
    v2df bcopy[nh][kk];
    v2df acopy[kk][kh];
    v2df ccopy[kk][kh];
    v2df acopy2[kk][kk];
    unsigned long bpcount, apcount, dotcount;
    bpcount= apcount= dotcount=0;
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy[j][k] = *((v2df*)(b[k]+j+j));
    //    printf("copy b end\n");
    for(i=0;i<m;i+=kk){
	for(k=0;k<kk;k++){
	    __builtin_prefetch(a+i+k+8,0,0);
	    __builtin_prefetch(c+i+k+8,1,0);
	    v2df  *ak = (v2df*)(a[i+k]);
	    v2df * awp = (v2df*)(acopy+k);
	    v2df  *ck = (v2df*)(c[i+k]);
	    v2df * cwp = (v2df*)(ccopy+k);
	    awp[0]=ak[0];
	    awp[1]=ak[1];
	    awp[2]=ak[2];
	    awp[3]=ak[3];
	    cwp[0]=ck[0];
	    cwp[1]=ck[1];
	    cwp[2]=ck[2];
	    cwp[3]=ck[3];
	}
	for (j=0;j<n;j++){
	    double * ap = (double*)( acopy+j);
	    for (k=0;k<kk;k++){
		acopy2[j][k]=(v2df){ap[k],ap[k]};
	    }
	}
	//	printf("copy a end\n");
	for(k=0;k<kk;k++){
	    v2df * cp = (v2df*) ccopy[k];
	    v2df * ap = acopy2[k];
	    for(j=0;j<nh;j++){
		v2df * bp = bcopy[j];
		cp[j] -= ap[0]*bp[0]
		    + ap[1]*bp[1]
		    + ap[2]*bp[2]
		    + ap[3]*bp[3]
		    + ap[4]*bp[4]
		    + ap[5]*bp[5]
		    + ap[6]*bp[6]
		    + ap[7]*bp[7];
	    }
	}
	for(k=0;k<kk;k++){
	    v2df  *ck = (v2df*)(c[i+k]);
	    v2df * cwp = (v2df*)(ccopy+k);
#if 0
	    ck[0] = cwp[0];
	    ck[1] = cwp[1];
	    ck[2] = cwp[2];
	    ck[3] = cwp[3];
#endif
	    __builtin_ia32_movntpd((double*)(ck),cwp[0]);
	    __builtin_ia32_movntpd((double*)(ck+1),cwp[1]);
	    __builtin_ia32_movntpd((double*)(ck+2),cwp[2]);
	    __builtin_ia32_movntpd((double*)(ck+3),cwp[3]);
	}
	//	printf("calc c end\n");
    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //	   (double)bpcount, (double)apcount, (double)dotcount);
}
void matmul_for_nk8(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{

    int i,j,k;
    const int kk = 8;
    const int kh = kk/2;
    int nh = n/2;
    v2df bcopy[nh][kk];
    v2df acopy2[kk][kk];
    //    unsigned long bpcount, apcount, dotcount;
    //    bpcount= apcount= dotcount=0;
    BEGIN_TSC;

    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy[j][k] = *((v2df*)(b[k]+j+j));
    //    END_TSC(bpcount);
    //    printf("copy b end\n");
#pragma omp parallel for private(i,j,k,acopy2)	 schedule(static)
    for(i=0;i<m;i+=kk){
	//	BEGIN_TSC;
	for(k=0;k<kk;k++){
	    __builtin_prefetch(a+i+k+16,0,0);
	    __builtin_prefetch(c+i+k+16,1,0);
	}
	for (j=0;j<n;j++){
	    double * ap = (double*)( a[i+j]);
	    for (k=0;k<kk;k++){
		acopy2[j][k]=(v2df){ap[k],ap[k]};
	    }
	}
	//	END_TSC(apcount);
	//	printf("copy a end\n");
	//	BEGIN_TSC;
	for(k=0;k<kk;k++){
	    v2df * cp = (v2df*) (c[i+k]);
	    v2df * ap = acopy2[k];
	    for(j=0;j<nh;j++){
		v2df * bp = bcopy[j];
		cp[j] -= ap[0]*bp[0] + ap[1]*bp[1]
		    + ap[2]*bp[2]    + ap[3]*bp[3]
		    + ap[4]*bp[4]    + ap[5]*bp[5]
		    + ap[6]*bp[6]    + ap[7]*bp[7];
	    }
	}
	//	printf("calc c end\n");
	//	END_TSC(dotcount);
    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //   (double)bpcount, (double)apcount, (double)dotcount);
    END_TSC(t,10);
}

void matmul_for_nk8_3(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{
    int ii;
    int dm = (m+31)/32;
    dm*= 8;
#pragma omp parallel for private(ii)	 schedule(static)
    for(ii=0;ii<4;ii++){
	int ifirst, iend;
	ifirst = ii*dm;
	iend = ifirst+dm;
	if (iend > m) iend = m;
	//	fprintf(stderr, "m, i, ifirst, iend = %d %d %d %d\n", m, ii, ifirst, iend);
	if (ifirst < m){
	    matmul_for_nk8_0(n1, a[ifirst], n2, b, n3, c[ifirst], iend-ifirst, n);
	}
    }

}

void matmul_for_nk16_0(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{
    int i;
    int mm=64;
    for(i=0;i<m;i+=mm){
	if (i+mm >m) mm=m-i;
	matmul_for_nk8_0(n1, (double(*)[]) (a[i]), n2, b,
			 n3, (double(*)[]) (c[i]), mm, 16);
	matmul_for_nk8_0(n1, (double(*)[]) (&a[i][8]), n2,(double(*)[])(b[8]),
			 n3, (double(*)[]) (c[i]), mm, 16);
    }

}
void matmul_for_nk16(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{
    if (m < 64){
	matmul_for_nk16_0c(n1, a, n2, b, n3, c, m, n);
	return;
    }
    int ii;
    int dm = (m+63)/64;
    dm*= 16;
#pragma omp parallel for private(ii)	 schedule(static)
    for(ii=0;ii<4;ii++){
	int ifirst, iend;
	ifirst = ii*dm;
	iend = ifirst+dm;
	if (iend > m) iend = m;
	//	fprintf(stderr, "m, i, ifirst, iend = %d %d %d %d\n", m, ii, ifirst, iend);
	if (ifirst < m){
	    matmul_for_nk16_0c(n1, a[ifirst], n2, b, n3, c[ifirst], iend-ifirst, n);
	}
    }

}

void matmul_for_nk32(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int n)
{
    int ii;
    int dm = (m+127)/128;
    dm*= 32;
#pragma omp parallel for private(ii)	 schedule(static)
    for(ii=0;ii<4;ii++){
	int ifirst, iend;
	ifirst = ii*dm;
	iend = ifirst+dm;
	if (iend > m) iend = m;
	//	fprintf(stderr, "m, i, ifirst, iend = %d %d %d %d\n", m, ii, ifirst, iend);
	if (ifirst < m){
	    matmul_for_nk32_0(n1, a[ifirst], n2, b, n3, c[ifirst], iend-ifirst, n);
	}
    }

}

void matmul_for_small_nk_7(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int kk,
			 int n)
{

    int i,j;
    int nh = n/2;
    register int k;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
    v2df acopy[kk];
    v2df acopy2[kk];
    double *acp = (double*) acopy;
    double *acp2 = (double*) acopy2;
    unsigned long bpcount, apcount, dotcount;
    if (kk == 8){
	matmul_for_nk8(n1,  a, n2, b, n3, c, m, n);
	return;
    }
    BEGIN_TSC;
    bpcount= apcount= dotcount=0;
    //    BEGIN_TSC;
    for(k=0;k<kk;k++)
	for(j=0;j<nh;j++)
		bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    //    END_TSC(bpcount);
    for(i=0;i<m;i+=2){
	//	BEGIN_TSC;
	double *ap=a[i];
	double *ap2=a[i+1];
	register v2df tmp, tmp2;
	v2df * cp = (v2df*) (&(c[i][0]));
	v2df * cp2 = (v2df*) (&(c[i+1][0]));

	for(k=0;k<kk;k+=2){
	    v2df * aa = (v2df*)(ap+k);
	     __builtin_prefetch((double*)a[i+4]+k,0);
	     __builtin_prefetch((double*)a[i+5]+k,0);
	    acopy [k]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	    acopy[k+1]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	    aa = (v2df*)(ap2+k);
	    acopy2 [k]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	    acopy2[k+1]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	    //   acp[k*2]=acp[k*2+1]=ap[k];
	    //	    acp[k*2+2]=acp[k*2+3]=ap[k+1];
	    //	    acp2[k*2]=acp2[k*2+1]=ap2[k];
	    //	    acp2[k*2+2]=acp2[k*2+3]=ap2[k+1];
	}
	//	END_TSC(apcount);
	//	BEGIN_TSC;
	for(j=0;j<nh;j++){
	    tmp = tmp2=  (v2df){0.0,0.0};
	    v2df ctmp= cp[j];
	    v2df ctmp2 = cp2[j] ;

	    v2df * bp = bcopy2[j];
	    __builtin_prefetch(c[i+4]+j,0);
	    __builtin_prefetch(c[i+5]+j,0);
	    for(k=0;k<kk;k+=8){
		int k2 = k+4;
		v2df *avp = acopy+k;
		v2df *avp2 = acopy2+k;
		v2df *bvp = bp+k;
		tmp += avp[0]*bvp[0];
		tmp2 += avp2[0]*bvp[0];
		tmp +=avp[1]*bvp[1];
		tmp2+=avp2[1]*bvp[1];
		tmp +=avp[2]*bvp[2];
		tmp2+=avp2[2]*bvp[2];
		tmp +=avp[3]*bvp[3];
		tmp2+=avp2[3]*bvp[3];
		tmp += avp[4]*bvp[4];
		tmp2 += avp2[4]*bvp[4];
		tmp +=avp[5]*bvp[5];
		tmp2+=avp2[5]*bvp[5];
		tmp +=avp[6]*bvp[6];
		tmp2+=avp2[6]*bvp[6];
		tmp +=avp[7]*bvp[7];
		tmp2+=avp2[7]*bvp[7];
	    }

	    cp[j] = ctmp -tmp;
	    cp2[j] = ctmp2 -tmp2;

	}
	//	END_TSC(dotcount);

    }
    //    printf("m, kk, n = %d %d %d counts = %g %g  %g\n", m,kk,n,
    //	   (double)bpcount, (double)apcount, (double)dotcount);
    END_TSC(t,11);

}

void matmul_for_small_nk(int n1, double a[][n1],
			 int n2, double b[][n2],
			 int n3, double c[][n3],
			 int m,
			 int kk,
			 int n)
{

    int i,ii;
    int nh = n/2;
    register int k;
    double bcopy[n][kk];
    v2df bcopy2[nh][kk];
    v2df acopy[kk];
    v2df acopy2[kk];
    double *acp = (double*) acopy;
    double *acp2 = (double*) acopy2;
    if (kk == 8){
	matmul_for_nk8_3(n1,  a, n2, b, n3, c, m, n);
	return;
    }
    if (kk == 16){
	matmul_for_nk16(n1,  a, n2, b, n3, c, m, n);
	return;
    }
    if (kk == 32){
	matmul_for_nk32(n1,  a, n2, b, n3, c, m, n);
	return;
    }
    BEGIN_TSC;
    for(k=0;k<kk;k++){
	int j;
	for(j=0;j<nh;j++)
	    bcopy2[j][k] = *((v2df*)(b[k]+j+j));
    }
#pragma omp parallel for private(i,k,acopy,acopy2)  schedule(static)
    for(i=0;i<m;i+=2){
	int j;
	double *ap=a[i];
	double *ap2=a[i+1];
	register v2df tmp, tmp2;
	v2df * cp = (v2df*) (&(c[i][0]));
	v2df * cp2 = (v2df*) (&(c[i+1][0]));
	for(k=0;k<kk;k+=2){
	    v2df * aa = (v2df*)(ap+k);
	    acopy [k]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	    acopy[k+1]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	    aa = (v2df*)(ap2+k);
	    acopy2 [k]= (v2df)  __builtin_ia32_shufpd(*aa,*aa,0x0);
	    acopy2[k+1]=(v2df) __builtin_ia32_shufpd(*aa,*aa,0xff);
	}
	__builtin_prefetch(a[i+4],0,3);
	__builtin_prefetch(c[i+4],1);
	__builtin_prefetch(a[i+5],0,3);
	__builtin_prefetch(c[i+5],1);
	__builtin_prefetch(a[i+20],0,3);
	__builtin_prefetch(c[i+20],1);
	__builtin_prefetch(a[i+21],0,3);
	__builtin_prefetch(c[i+21],1);
	for(j=0;j<nh;j++){
	    tmp = tmp2=  (v2df){0.0,0.0};
	    v2df ctmp= cp[j];
	    v2df ctmp2 = cp2[j] ;
	    v2df * bp = bcopy2[j];
	    for(k=0;k<kk;k+=8){
		int k2 = k+4;
		v2df *avp = acopy+k;
		v2df *avp2 = acopy2+k;
		v2df *bvp = bp+k;
		tmp += avp[0]*bvp[0];
		tmp2 += avp2[0]*bvp[0];
		tmp +=avp[1]*bvp[1];
		tmp2+=avp2[1]*bvp[1];
		tmp +=avp[2]*bvp[2];
		tmp2+=avp2[2]*bvp[2];
		tmp +=avp[3]*bvp[3];
		tmp2+=avp2[3]*bvp[3];
		tmp += avp[4]*bvp[4];
		tmp2 += avp2[4]*bvp[4];
		tmp +=avp[5]*bvp[5];
		tmp2+=avp2[5]*bvp[5];
		tmp +=avp[6]*bvp[6];
		tmp2+=avp2[6]*bvp[6];
		tmp +=avp[7]*bvp[7];
		tmp2+=avp2[7]*bvp[7];
	    }
	    cp[j] = ctmp -tmp;
	    cp2[j] = ctmp2 -tmp2;
	}
    }
    END_TSC(t,11);
}

void mydgemm(int m,
	     int n,
	     int k,
	     double alpha,
	     double * a,
	     int na,
	     double * b,
	     int nb,
	     double beta,
	     double * c,
	     int nc)
{
    double t0, t1, t2;
    //    printf("mydgemm called %d %d %d\n", m, n, k);
    if (k>= 512){
	get_cputime(&t0,&t1);
    }
    BEGIN_TSC;
    BEGIN_TIMER(timer);
#ifdef USEGDR
    if ((k>512)  || ((k==512) && ((n>=1024)||(m>=1024)))){
	//    if (k>=2048){
	mygdrdgemm(m, n, k, alpha, a, na, b, nb, beta, c, nc);
    }else{
	if ((k<=16) && (alpha == -1.0) && (beta == 1.0)){
	    matmul_for_small_nk(na, a,  nb, b, nc, c,  m, n, k);
	}else{

	    cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans,
			 m,n, k, alpha, a, na, b, nb, beta, c, nc);
	}
    }
#else
#ifdef GEMMTEST
    mytestdgemm(m, n, k, na, nb, nc,alpha, a, b, beta, c);
#else
    cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans,
	         m,n, k, alpha, a, na, b, nb, beta, c, nc);
#endif
#endif
    //    cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans,
    //		 m,n, k, alpha, a, na, b, nb, beta, c, nc);
    if (k==2048){
	END_TIMER(timer,31,((double)(m))*n*k*2);
	END_TSC(t,14);
    }else if (k==1024){
	END_TIMER(timer,32,((double)(m))*n*k*2);
	END_TSC(t,15);
    }else if (k==512){
	END_TIMER(timer,33,((double)(m))*n*k*2);
	END_TSC(t,17);
    }else{
	END_TIMER(timer,34,((double)(m))*n*k*2);
	END_TSC(t,18);
    }

    if (k>= 512){
	get_cputime(&t0,&t1);
	dprintf(10,"dgemm M=%d N=%d K=%d time=%10.4g  %g Gflops\n",
		m,n,k,t0, ((double)m)*n*k*2/t0/1e9);
    }
}

void reset_gdr(int m, double a[][m], int nb, double awork[][nb], int n)
{
#ifdef USEGDR
    double aw2[nb][nb];
    if (nb < 2048){
	fprintf(stderr,"reset_gdr nb = %d <2048 not supported\n", nb);
	exit(-1);
    }
    gdr_check_and_restart(a, awork, aw2);

    int i,j;
    dprintf(9,"reset_gdr clear awork\n");
    for (i=0;i<nb;i++){
	for (j=0;j<n;j++){
	    awork[j][i]=0;
	}
    }
    dprintf(9,"reset_gdr clear aw2\n");
    for (i=0;i<nb;i++){
	for (j=0;j<nb;j++){
	    aw2[j][i]=0;
	}
    }
    dprintf(9,"reset_gdr try_dgemm\n");
    gdrsetforceswapab();
    cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans,
	         n,nb, nb, 1.0, awork, nb, aw2, nb, 0.0, a, m);
    mydgemm(n,nb,nb,1.0,awork,nb,aw2,nb,0.0,a,m);

#endif
}

#ifndef USEGDR
void gdrsetforceswapab(){}
void gdrresetforceswapab(){}
void gdrsetskipsendjmat(){};
void gdrresetskipsendjmat(){}
void gdrsetnboards(){}
void set_matmul_msg_level(int level){}
void gdrdgemm_set_stress_factor(int x){}
#endif


text:
/*
#
# Copyright (c) 2006-2012      University of Houston. All rights reserved.
# $COPYRIGHT$
#
# Additional copyrights may follow
#
# $HEADER$
#
*/

#ifndef __SL_INTERNAL__
#define __SL_INTERNAL__

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <sys/types.h>
//#include <gsl/gsl_fit.h>
#include <sys/time.h>
#include <time.h>
#include <unistd.h>
#include <fcntl.h>
#include <math.h>
#include <errno.h>
//#include <papi.h>
#ifdef MINGW
#include <windows.h>
#include <winsock2.h>
#include <ws2tcpip.h>
#else
#include <sys/socket.h>
#include <netinet/in.h>
#include <netinet/tcp.h>
#include <arpa/inet.h>
#include <netdb.h>
#include <sys/wait.h>
#include <sys/resource.h>
#include <pwd.h>
#include <sys/utsname.h>
#endif

#include "SL_array.h"

extern fd_set SL_send_fdset;
extern fd_set SL_recv_fdset;

extern int SL_this_procid;
extern int SL_this_procport;
extern int SL_this_listensock;
extern int SL_numprocs;
extern int SL_init_numprocs;
extern int SL_proxy_numprocs;
/* Message header send before any message */
struct SL_msg_header {
    int      cmd; /* what type of message is this */
    int     from; /* id of the src process */
    int       to; /* id of the dest. process */
    int      tag; /* tag of the message */
    int  context; /* context id */
    int      len; /* Message length in bytes */
    int       id; /* Id of the last fragment */
    int loglength;
    int     temp;
};
typedef struct SL_msg_header SL_msg_header;

/* Process structure containing all relevant contact information,
   communication status etc. */
struct SL_proc;
typedef int SL_msg_comm_fnct ( struct SL_proc *dproc, int fd );

struct SL_proc {
    int                         id;
    char                 *hostname;
    int                       port;
    int                       sock;
    int                      state;
    int           connect_attempts; /* number of connect attempts */
    double    connect_start_tstamp; /* time stamp when we started to accept or connect
				       for this proc */
    double                 timeout; /* max time a process should wait before disconnecting */

    struct SL_msgq_head    *squeue; /* Send queue */
    struct SL_msgq_head    *rqueue; /* Recv queue */
    struct SL_msgq_head   *urqueue; /* Unexpected msgs queue */
    struct SL_msgq_head   *scqueue; /* Send complete queue */
    struct SL_msgq_head   *rcqueue; /* Recv complete queue */
    struct SL_qitem   *currecvelem;
    struct SL_qitem   *cursendelem;
    SL_msg_comm_fnct     *recvfunc;
    SL_msg_comm_fnct     *sendfunc;
    struct SL_msg_perf    *msgperf; /*to keep track of time and msglenth for each communication */
    struct SL_msg_perf   *insertpt;
    struct SL_network_perf *netperf;
};

struct SL_msg_perf {
	struct SL_msg_perf  *fwd;
	struct SL_msg_perf  *back;
	int 		   msglen;
	double		     time;
	int 		      pos;
	int 		  msgtype; /*send type(0) or recieve type(1)*/
	int 		   elemid;
        struct SL_proc      *proc;
};
typedef struct SL_msg_perf SL_msg_perf;

struct SL_network_perf {
	struct SL_network_perf	 *fwd;
	struct SL_network_perf	*back;
	double 	              latency;
	double 		    bandwidth;
	int 			  pos;

};
typedef struct SL_network_perf SL_network_perf;

typedef struct SL_proc SL_proc;

#ifdef MINGW
struct iovec {
	char	*iov_base;	/* Base address. */
	size_t	 iov_len;	/* Length. */
};

#define	TCP_MAXSEG		0x02	/* set maximum segment size */
#define 	F_GETFL		3	/* get file->f_flags */
#define 	F_SETFL		4	/* set file->f_flags */
#define 	O_NONBLOCK	 	00004
#endif

/* A  message queue item containing the operation
   it decsribes */
struct SL_msgq_head;
struct SL_qitem {
    int                       id;
    int                   iovpos;
    int                   lenpos;
    int                    error;
    struct iovec          iov[2];
    struct SL_msgq_head *move_to;
    struct SL_msgq_head    *head;
    struct SL_qitem        *next;
    struct SL_qitem        *prev;
    double             starttime;
    double               endtime;
};
typedef struct SL_qitem SL_qitem;

/* A message queue */
struct SL_msgq_head {
    int               count;
    char              *name;
    struct SL_qitem  *first;
    struct SL_qitem   *last;
};
typedef struct SL_msgq_head SL_msgq_head;

/* Request object identifying an ongoing communication */
struct SL_msg_request {
    struct SL_proc        *proc;
    int                    type; /* Send or Recv */
    int                      id;
    struct SL_qitem       *elem;
    struct SL_msgq_head *cqueue; /* completion queue to look for */
};
typedef struct SL_msg_request SL_msg_request;

struct SL_msgq_head     *SL_event_sendq;
struct SL_msgq_head     *SL_event_recvq;
struct SL_msgq_head     *SL_event_sendcq;

/* MACROS */
/*#ifdef PRINTF
  #undef PRINTF
  #define PRINTF(A) printf A
#else
  #define PRINTF(A)
#endif*/

#define FALSE 0
#define TRUE  1

#define SEND 0
#define RECV 1

#define SL_RECONN_MAX      20
#define SL_ACCEPT_MAX_TIME 10
#define SL_READ_MAX_TIME   5
#define SL_ACCEPT_INFINITE_TIME -1
#define SL_BIND_PORTSHIFT 200
#define SL_SLEEP_TIME       1
#define SL_TCP_BUFFER_SIZE  262142
#define SL_MAX_EVENT_HANDLE     1
#define SL_CONSTANT_ID          -32
#define SL_EVENT_MANAGER        -1
#define SL_PROXY_SERVER		-2
#define PERFBUFSIZE	20
#define MTU		(1024L*4L)
#define SL_PROC_ID	-64

int SL_socket ( void );
int SL_bind_static ( int handle, int port );
int SL_bind_dynamic (  int handle, int port );
int SL_socket_close ( int handle );
int SL_open_socket_conn ( int *handle, const char *as_host, int port );
int SL_open_socket_bind ( int *handle, int port );
int SL_open_socket_listen ( int sock );
int SL_open_socket_listen_nb ( int *handle, int port );
int SL_open_socket_conn_nb ( int *handle, const char *as_host, int port );
int SL_socket_read ( int hdl, char *buf, int num, double timeout );
int SL_socket_write ( int hdl, char *buf, int num, double timeout );
int SL_socket_write_nb ( int hdl, char *buf, int num, int *numwritten );

int SL_socket_read_nb ( int hdl, char *buf, int num, int* numread );

void SL_print_socket_options ( int fd );
void SL_configure_socket ( int sd );
void SL_configure_socket_nb ( int sd );
int SL_init_internal();
double SL_papi_time();

/* status object t.b.d */

#endif /* __SL_INTERNALL__ */

ext: h

max_stars_repo_name: edgargabriel/VolpexMPI

max_issues_repo_name: edgargabriel/VolpexMPI

max_forks_repo_name: edgargabriel/VolpexMPI


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef delaunay_3673850c_175d_4077_bef2_fc69e19e5cf8_h
#define delaunay_3673850c_175d_4077_bef2_fc69e19e5cf8_h

#include <gslib/std.h>
#include <gslib/math.h>
#include <ariel/config.h>

/*
 * A implementation of the Delaunay Triangulation and the Constrained Delaunay Triangulation.
 * The original implementation was:
 * https://code.google.com/archive/p/cgprojects09/
 * further more, modifications were done to achieve the constrained version, visit:
 * http://www.geom.uiuc.edu/~samuelp/del_project.html
 */

__ariel_begin__

class ariel_export dt_joint
{
protected:
    vec2            _point;
    void*           _binding;

public:
    dt_joint() { _binding = nullptr; }
    dt_joint(const vec2& p, void* b) { _point = p, _binding = b; }
    void set_point(const vec2& p) { _point = p; }
    const vec2& get_point() const { return _point; }
    void set_binding(void* b) { _binding = b; }
    void* get_binding() const { return _binding; }
};

class ariel_export dt_edge
{
protected:
    dt_joint*       _org;
    dt_edge*        _prev;
    dt_edge*        _next;
    dt_edge*        _symmetric;
    bool            _constraint;
    bool            _checked;

public:
    dt_edge();
    void set_org(dt_joint* j) { _org = j; }
    dt_joint* get_org() const { return _org; }
    void set_dest(dt_joint* j) { _symmetric->set_org(j); }
    dt_joint* get_dest() const { return _symmetric->get_org(); }
    void set_prev_edge(dt_edge* e) { _prev = e; }
    dt_edge* get_prev_edge() const { return _prev; }
    void set_next_edge(dt_edge* e) { _next = e; }
    dt_edge* get_next_edge() const { return _next; }
    void set_symmetric(dt_edge* e) { _symmetric = e; }
    dt_edge* get_symmetric() const { return _symmetric; }
    const vec2& get_org_point() const { return _org->get_point(); }
    const vec2& get_dest_point() const { return get_dest()->get_point(); }
    void* get_org_binding() const { return _org->get_binding(); }
    void* get_dest_binding() const { return get_dest()->get_binding(); }
    dt_edge* get_org_next() const { return _prev->get_symmetric(); }
    dt_edge* get_org_prev() const { return _symmetric->get_next_edge(); }
    dt_edge* get_dest_next() const { return _symmetric->get_prev_edge(); }
    dt_edge* get_dest_prev() const { return _next->get_symmetric(); }
    dt_edge* get_left_next() const { return _next; }
    dt_edge* get_left_prev() const { return _prev; }
    dt_edge* get_right_next() const { return _symmetric->get_next_edge()->get_symmetric(); }
    dt_edge* get_right_prev() const { return _symmetric->get_prev_edge()->get_symmetric(); }
    void set_constraint(bool b) { _constraint = b; }
    bool is_constraint() const { return _constraint; }
    void set_checked(bool c) { _checked = c; }
    bool is_checked() const { return _checked; }
    bool is_outside_boundary() const;
    bool is_boundary() const;
};

struct ariel_export dt_edge_range
{
    dt_edge*        left;
    dt_edge*        right;

    dt_edge_range()  { left = right = nullptr; }
};

typedef list<dt_joint> dt_input_joints;
typedef list<dt_joint*> dt_joint_ptrs;
typedef vector<dt_joint*> dt_joints;
typedef unordered_set<dt_edge*> dt_edges;
typedef vector<dt_edge*> dt_edge_list;

class ariel_export delaunay_triangulation
{
public:
    delaunay_triangulation() {}
    ~delaunay_triangulation() { clear(); }
    void initialize(dt_input_joints& inputs);
    void run();
    void clear();
    dt_edge* add_constraint(const vec2& p1, const vec2& p2);
    void trim(dt_edge_list& edges);
    void set_range_left(dt_edge* e) { _edge_range.left = e; }
    void tracing() const;
    void trace_heuristically() const;
    void trace_mel() const;

protected:
    dt_joints       _sorted_joints;
    dt_edge_range   _edge_range;
    dt_edges        _edge_holdings;

protected:
    bool is_in_range(int begin, int end, dt_joint* joint);
    dt_edge_range delaunay(int begin, int end);
    dt_edge* create_edge_pair();
    dt_edge* connect_edges(dt_edge* e1, dt_edge* e2);
    void destroy_edge_pair(dt_edge* e);
    void shrink_triangulate(dt_edge* cut);
    void shrink_recursively(dt_edge_list& strips, dt_edge_list& temps);
    void collect_trim_edges(dt_edge_list& edges, dt_edge* e);

    template<class _visit>
    static void traverse_triangles_recursively(dt_edge* e, _visit visit)
    {
        assert(e);
        if(e->is_checked())
            return;
        auto do_visit = [&visit](dt_edge* e) {
            assert(e);
            auto* e1 = e->get_prev_edge();
            auto* e2 = e->get_next_edge();
            assert(e1 && e2 && !e1->is_checked() && !e2->is_checked());
            void* b1 = e->get_org_binding();
            void* b2 = e->get_dest_binding();
            void* b3 = e1->get_org_binding();
            void* b4 = e2->get_dest_binding();
            if(b3 != b4)    /* invalid triangle */
                return;
            auto& p1 = e->get_org_point();
            auto& p2 = e->get_dest_point();
            auto& p3 = e1->get_org_point();
            if(!is_concave_angle(p1, p2, p3))     /* ccw! */
                return;
            bool b[3];
            b[0] = e->is_boundary();
            b[1] = e1->is_boundary();
            b[2] = e2->is_boundary();
            visit(b1, b2, b3, b);
            e->set_checked(true);
            e1->set_checked(true);
            e2->set_checked(true);
        };
        do_visit(e);
        traverse_triangles_recursively(e->get_symmetric(), visit);
        traverse_triangles_recursively(e->get_prev_edge()->get_symmetric(), visit);
        traverse_triangles_recursively(e->get_next_edge()->get_symmetric(), visit);
    }
    template<class _visit>
    static void traverse_triangles_flatly(dt_edge* e, _visit visit)
    {
        vector<dt_edge*> edges;
        edges.push_back(e);
        while(!edges.empty()) {
            e = edges.back();
            edges.pop_back();
            if(e->is_checked())
                continue;
            auto do_visit = [&visit](dt_edge* e) {
                assert(e);
                auto* e1 = e->get_prev_edge();
                auto* e2 = e->get_next_edge();
                assert(e1 && e2 && !e1->is_checked() && !e2->is_checked());
                void* b1 = e->get_org_binding();
                void* b2 = e->get_dest_binding();
                void* b3 = e1->get_org_binding();
                void* b4 = e2->get_dest_binding();
                if(b3 != b4)    /* invalid triangle */
                    return;
                auto& p1 = e->get_org_point();
                auto& p2 = e->get_dest_point();
                auto& p3 = e1->get_org_point();
                if(!is_concave_angle(p1, p2, p3))     /* ccw! */
                    return;
                bool b[3];
                b[0] = e->is_boundary();
                b[1] = e1->is_boundary();
                b[2] = e2->is_boundary();
                visit(b1, b2, b3, b);
                e->set_checked(true);
                e1->set_checked(true);
                e2->set_checked(true);
            };
            do_visit(e);
            edges.push_back(e->get_symmetric());
            edges.push_back(e->get_prev_edge()->get_symmetric());
            edges.push_back(e->get_next_edge()->get_symmetric());
        }
    }

public:
    template<class _visit>
    void traverse_triangles(_visit visit)
    {
        traverse_triangles_flatly(_edge_range.left, visit);
    }
    void reset_traverse();
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/* specfunc/beta_inc.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 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, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */
/* Modified for cdfs by Brian Gough, June 2003 */

#include <gsl/gsl_sf_gamma.h>

static double
beta_cont_frac (const double a, const double b, const double x,
                const double epsabs)
{
  const unsigned int max_iter = 512;    /* control iterations      */
  const double cutoff = 2.0 * GSL_DBL_MIN;      /* control the zero cutoff */
  unsigned int iter_count = 0;
  double cf;

  /* standard initialization for continued fraction */
  double num_term = 1.0;
  double den_term = 1.0 - (a + b) * x / (a + 1.0);

  if (fabs (den_term) < cutoff)
    den_term = GSL_NAN;

  den_term = 1.0 / den_term;
  cf = den_term;

  while (iter_count < max_iter)
    {
      const int k = iter_count + 1;
      double coeff = k * (b - k) * x / (((a - 1.0) + 2 * k) * (a + 2 * k));
      double delta_frac;

      /* first step */
      den_term = 1.0 + coeff * den_term;
      num_term = 1.0 + coeff / num_term;

      if (fabs (den_term) < cutoff)
        den_term = GSL_NAN;

      if (fabs (num_term) < cutoff)
        num_term = GSL_NAN;

      den_term = 1.0 / den_term;

      delta_frac = den_term * num_term;
      cf *= delta_frac;

      coeff = -(a + k) * (a + b + k) * x / ((a + 2 * k) * (a + 2 * k + 1.0));

      /* second step */
      den_term = 1.0 + coeff * den_term;
      num_term = 1.0 + coeff / num_term;

      if (fabs (den_term) < cutoff)
        den_term = GSL_NAN;

      if (fabs (num_term) < cutoff)
        num_term = GSL_NAN;

      den_term = 1.0 / den_term;

      delta_frac = den_term * num_term;
      cf *= delta_frac;

      if (fabs (delta_frac - 1.0) < 2.0 * GSL_DBL_EPSILON)
        break;

      if (cf * fabs (delta_frac - 1.0) < epsabs)
        break;

      ++iter_count;
    }

  if (iter_count >= max_iter)
    return GSL_NAN;

  return cf;
}

/* The function beta_inc_AXPY(A,Y,a,b,x) computes A * beta_inc(a,b,x)
   + Y taking account of possible cancellations when using the
   hypergeometric transformation beta_inc(a,b,x)=1-beta_inc(b,a,1-x).

   It also adjusts the accuracy of beta_inc() to fit the overall
   absolute error when A*beta_inc is added to Y. (e.g. if Y >>
   A*beta_inc then the accuracy of beta_inc can be reduced) */

static double
beta_inc_AXPY (const double A, const double Y,
               const double a, const double b, const double x)
{
  if (x == 0.0)
    {
      return A * 0 + Y;
    }
  else if (x == 1.0)
    {
      return A * 1 + Y;
    }
  else if (a > 1e5 && b < 10 && x > a / (a + b))
    {
      /* Handle asymptotic regime, large a, small b, x > peak [AS 26.5.17] */
      double N = a + (b - 1.0) / 2.0;
      return A * gsl_sf_gamma_inc_Q (b, -N * log (x)) + Y;
    }
  else if (b > 1e5 && a < 10 && x < b / (a + b))
    {
      /* Handle asymptotic regime, small a, large b, x < peak [AS 26.5.17] */
      double N = b + (a - 1.0) / 2.0;
      return A * gsl_sf_gamma_inc_P (a, -N * log1p (-x)) + Y;
    }
  else
    {
      double ln_beta = gsl_sf_lnbeta (a, b);
      double ln_pre = -ln_beta + a * log (x) + b * log1p (-x);

      double prefactor = exp (ln_pre);

      if (x < (a + 1.0) / (a + b + 2.0))
        {
          /* Apply continued fraction directly. */
          double epsabs = fabs (Y / (A * prefactor / a)) * GSL_DBL_EPSILON;

          double cf = beta_cont_frac (a, b, x, epsabs);

          return A * (prefactor * cf / a) + Y;
        }
      else
        {
          /* Apply continued fraction after hypergeometric transformation. */
          double epsabs =
            fabs ((A + Y) / (A * prefactor / b)) * GSL_DBL_EPSILON;
          double cf = beta_cont_frac (b, a, 1.0 - x, epsabs);
          double term = prefactor * cf / b;

          if (A == -Y)
            {
              return -A * term;
            }
          else
            {
              return A * (1 - term) + Y;
            }
        }
    }
}

/* Direct series evaluation for testing purposes only */

#if 0
static double
beta_series (const double a, const double b, const double x,
             const double epsabs)
{
  double f = x / (1 - x);
  double c = (b - 1) / (a + 1) * f;
  double s = 1;
  double n = 0;

  s += c;

  do
    {
      n++;
      c *= -f * (2 + n - b) / (2 + n + a);
      s += c;
    }
  while (n < 512 && fabs (c) > GSL_DBL_EPSILON * fabs (s) + epsabs);

  s /= (1 - x);

  return s;
}
#endif

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_interp.h>

typedef double (*eval_func)(const gsl_interp*,const double *x,const double *y,
                            double,gsl_interp_accel*);

static const double EPSILON = 0.00000001;
static const double EPS = 0.0001;
static const int PTS = 801;

static const double x_values[] = { 0, 1, 2, 3, 4, 5 };
static const double y_values[] = { 5, 0, 5, 2, 0, 5 };
static double x_range;

static const int N = sizeof(x_values)/sizeof(double);

static double* interpolate(const gsl_interp_type *type, const size_t order);
static void    check_discont(gsl_interp *interp, const eval_func func,
                             const char *dsym, const double x,
                             gsl_interp_accel *accel);

static inline double
scale_x(const int i)
{
  return x_values[0] + (double)i/(PTS-1)*x_range;
}

int
main(int argc, char *argv[])
{
  int i;
  double *cspline, *cspline_d, *cspline_dd;
  double *akima, *akima_d, *akima_dd;

  x_range = x_values[N-1] - x_values[0];

  akima = interpolate(gsl_interp_akima_periodic, 0);
  akima_d = interpolate(gsl_interp_akima_periodic, 1);
  akima_dd = interpolate(gsl_interp_akima_periodic, 2);
  cspline = interpolate(gsl_interp_cspline_periodic, 0);
  cspline_d = interpolate(gsl_interp_cspline_periodic, 1);
  cspline_dd = interpolate(gsl_interp_cspline_periodic, 2);

  puts("# x   akima  cspline  akima'  cspline'  akima''  cspline''");
  for (i = -PTS/2; i < PTS/2; i++)
    printf("%0.5f %0.5f %0.5f %0.5f %0.5f %0.5f %0.5f\n",
           scale_x(i),
           akima[(i + PTS)%PTS],
           cspline[(i + PTS)%PTS],
           akima_d[(i + PTS)%PTS],
           cspline_d[(i + PTS)%PTS],
           akima_dd[(i + PTS)%PTS],
           cspline_dd[(i + PTS)%PTS]);

  free(akima);
  free(akima_d);
  free(akima_dd);
  free(cspline);
  free(cspline_d);
  free(cspline_dd);
  return 0;
}

static double*
interpolate(const gsl_interp_type *type, const size_t order)
{
  static const eval_func DERIV_ORDERS[] = {
    gsl_interp_eval,
    gsl_interp_eval_deriv,
    gsl_interp_eval_deriv2
  };
  static const char *DERIV_SYMS[] = {
    "",
    "'",
    "''"
  };

  gsl_interp *interp = gsl_interp_alloc(type, N);
  gsl_interp_accel *accel = gsl_interp_accel_alloc();
  double *plotpts = (double*)malloc(PTS*sizeof(double));
  int i;

  assert(order < sizeof(DERIV_ORDERS));
  gsl_interp_init(interp, x_values, y_values, N);
  for (i = 0; i < PTS; i++)
    plotpts[i] = DERIV_ORDERS[order](interp, x_values, y_values,
                                     scale_x(i), accel);

  for (i = 0; i < N-1; i++)
    check_discont(interp, DERIV_ORDERS[order], DERIV_SYMS[order],
                  x_values[i], accel);

  gsl_interp_free(interp);
  gsl_interp_accel_free(accel);

  return plotpts;
}

static void
check_discont(gsl_interp *interp, const eval_func func, const char *dsym,
              const double x, gsl_interp_accel *accel)
{
  double xm, xp;
  double ym, yp;

  xm = fmod(x - EPSILON + x_range, x_range);
  xp = fmod(x + EPSILON, x_range);
  ym = func(interp, x_values, y_values, xm, accel);
  yp = func(interp, x_values, y_values, xp, accel);
  if (fabs(yp - ym) > EPS) {
    fprintf(stderr, "discontinuity in %s%s at %f:\n"
                    "\t left value: %f\n"
                    "\tright value: %f\n"
                    "\t difference: %f\n\n",
            gsl_interp_name(interp), dsym, x, ym, yp, yp-ym);
  }
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
//
// Created by Harold on 2020/9/16.
//

#ifndef M_MATH_M_INTERP_H
#define M_MATH_M_INTERP_H

#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include <vector>
#include <algorithm>

namespace M_MATH {
    class Interpolation1D {
    public:
        enum InterpType {
            LINEAR,
            POLYNOMIAL,
            CSPLINE,
            CSPLINE_PERIODIC,
            AKIMA,
            AKIMA_PERIODIC,
            STEFFEN
        };

        static bool interpolate(double const* x,
                                double const* y,
                                size_t length,
                                InterpType type,
                                double const* n_x,
                                size_t n_length,
                                double* n_y);

        // this will create tmp var inside, better to use former one if possible
        template<typename IT_X, typename IT_Y, typename IT_NX, typename IT_NY>
        static bool interpolate(IT_X x_begin, IT_X x_end,
                                IT_Y y_begin, IT_Y y_end,
                                InterpType type,
                                IT_NX nx_begin, IT_NX nx_end,
                                IT_NY ny_begin);

    private:
        static const gsl_interp_type* InterpType2GSL(InterpType type);
    };

    class Interpolation2D {
    public:
        enum InterpType {
            BI_LINEAR,
            BI_CUBIC
        };

        static bool interpolate(double const* x,
                                double const* y,
                                double const* z,
                                size_t x_length,
                                size_t y_length,
                                InterpType type,
                                double const* n_x,
                                double const* n_y,
                                size_t n_x_length,
                                size_t n_y_length,
                                double* n_z);

        // this will create tmp var inside, better to use former one if possible
        template<typename IT_X, typename IT_Y, typename IT_Z, typename IT_NX, typename IT_NY, typename IT_NZ>
        static bool interpolate(IT_X x_begin, IT_X x_end,
                                IT_Y y_begin, IT_Y y_end,
                                IT_Z z_begin, IT_Z z_end,
                                InterpType type,
                                IT_NX nx_begin, IT_NX nx_end,
                                IT_NY ny_begin, IT_NY ny_end,
                                IT_NZ nz_begin);

    private:
        static const gsl_interp2d_type* InterpType2GSL(InterpType type);
    };

    const gsl_interp_type *Interpolation1D::InterpType2GSL(Interpolation1D::InterpType type) {
        switch (type) {
            case LINEAR:
                return gsl_interp_linear;
            case POLYNOMIAL:
                return gsl_interp_polynomial;
            case CSPLINE:
                return gsl_interp_cspline;
            case CSPLINE_PERIODIC:
                return gsl_interp_cspline_periodic;
            case AKIMA:
                return gsl_interp_akima;
            case AKIMA_PERIODIC:
                return gsl_interp_akima_periodic;
            case STEFFEN:
                return gsl_interp_steffen;
        }
        return nullptr;
    }

    bool Interpolation1D::interpolate(const double *x,
                                      const double *y,
                                      size_t length,
                                      Interpolation1D::InterpType type,
                                      const double *n_x,
                                      size_t n_length,
                                      double* n_y) {
        gsl_spline* spl = gsl_spline_alloc(InterpType2GSL(type), length);
        if (spl == nullptr)
            return false;
        gsl_interp_accel* acc = gsl_interp_accel_alloc();
        if (acc == nullptr) {
            return false;
        }
        gsl_spline_init(spl, x, y, length);
        for (auto i = 0; i < n_length; ++i)
            n_y[i] = gsl_spline_eval(spl, n_x[i], acc);
        gsl_spline_free(spl);
        gsl_interp_accel_free(acc);
        return true;
    }

    template<typename IT_X, typename IT_Y, typename IT_NX, typename IT_NY>
    bool Interpolation1D::interpolate(const IT_X x_begin, const IT_X x_end,
                                      const IT_Y y_begin, const IT_Y y_end,
                                      Interpolation1D::InterpType type,
                                      const IT_NX nx_begin, const IT_NX nx_end,
                                      IT_NY ny_begin) {
        auto len = std::min(std::distance(x_begin, x_end), std::distance(y_begin, y_end));
        auto n_len = std::distance(nx_begin, nx_end);
        std::vector<double> x_tmp(len), y_tmp(len);
        std::copy(x_begin, x_begin+len, x_tmp.begin());
        std::copy(y_begin, y_begin+len, y_tmp.begin());

        gsl_spline* spl = gsl_spline_alloc(InterpType2GSL(type), len);
        if (spl == nullptr)
            return false;
        gsl_interp_accel* acc = gsl_interp_accel_alloc();
        if (acc == nullptr) {
            return false;
        }
        gsl_spline_init(spl, x_tmp.data(), y_tmp.data(), len);
        for (auto i = 0; i < n_len; ++i)
            *(ny_begin+i) = gsl_spline_eval(spl, *(nx_begin+i), acc);
        gsl_spline_free(spl);
        gsl_interp_accel_free(acc);
        return true;
    }

    const gsl_interp2d_type *Interpolation2D::InterpType2GSL(Interpolation2D::InterpType type) {
        switch (type) {
            case BI_LINEAR:
                return gsl_interp2d_bilinear;
            case BI_CUBIC:
                return gsl_interp2d_bicubic;
        }
        return nullptr;
    }

    bool Interpolation2D::interpolate(const double *x,
                                      const double *y,
                                      const double *z,
                                      size_t x_length,
                                      size_t y_length,
                                      Interpolation2D::InterpType type,
                                      const double *n_x,
                                      const double *n_y,
                                      size_t n_x_length,
                                      size_t n_y_length,
                                      double *n_z) {
        gsl_spline2d* spl = gsl_spline2d_alloc(InterpType2GSL(type), x_length, y_length);
        if (spl == nullptr)
            return false;
        gsl_interp_accel *x_acc = gsl_interp_accel_alloc();
        if (x_acc == nullptr)
            return false;
        gsl_interp_accel *y_acc = gsl_interp_accel_alloc();
        if (y_acc == nullptr)
            return false;
        // zij = z[j * x_length + i]
        gsl_spline2d_init(spl, x, y, z, x_length, y_length);
        for (auto i = 0; i < n_x_length; ++i)
            for (auto j = 0; j < n_y_length; ++j)
                n_z[j * n_x_length + i] = gsl_spline2d_eval(spl, n_x[i], n_y[j], x_acc, y_acc);
        gsl_spline2d_free(spl);
        gsl_interp_accel_free(x_acc);
        gsl_interp_accel_free(y_acc);
        return true;
    }

    template<typename IT_X, typename IT_Y, typename IT_Z, typename IT_NX, typename IT_NY, typename IT_NZ>
    bool Interpolation2D::interpolate(const IT_X x_begin, const IT_X x_end,
                                      const IT_Y y_begin, const IT_Y y_end,
                                      const IT_Z z_begin, const IT_Z z_end,
                                      Interpolation2D::InterpType type,
                                      const IT_NX nx_begin, const IT_NX nx_end,
                                      const IT_NY ny_begin, const IT_NY ny_end,
                                      const IT_NZ nz_begin) {
        auto x_length = std::distance(x_begin, x_end);
        auto y_length = std::distance(y_begin, y_end);
        auto z_length = std::distance(z_begin, z_end);
        auto n_x_length = std::distance(nx_begin, nx_end);
        auto n_y_length = std::distance(ny_begin, ny_end);
        std::vector<double> x_tmp(x_length), y_tmp(y_length), z_tmp(z_length);
        std::copy(x_begin, x_end, x_tmp.begin());
        std::copy(y_begin, y_end, y_tmp.begin());
        std::copy(z_begin, z_end, z_tmp.begin());

        gsl_spline2d* spl = gsl_spline2d_alloc(InterpType2GSL(type), x_length, y_length);
        if (spl == nullptr)
            return false;
        gsl_interp_accel *x_acc = gsl_interp_accel_alloc();
        if (x_acc == nullptr)
            return false;
        gsl_interp_accel *y_acc = gsl_interp_accel_alloc();
        if (y_acc == nullptr)
            return false;
        // zij = z[j * x_length + i]
        gsl_spline2d_init(spl, x_tmp.data(), y_tmp.data(), z_tmp.data(), x_length, y_length);
        for (auto i = 0; i < n_x_length; ++i)
            for (auto j = 0; j < n_y_length; ++j)
                *(nz_begin + j * n_x_length + i) = gsl_spline2d_eval(spl, *(nx_begin + i), *(ny_begin + j), x_acc, y_acc);
        gsl_spline2d_free(spl);
        gsl_interp_accel_free(x_acc);
        gsl_interp_accel_free(y_acc);
        return true;
    }
}

#endif //M_MATH_M_INTERP_H


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <fcntl.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <sys/file.h>
#include <unistd.h>
#include <gsl/gsl_rng.h>

#include "allvars.h"
#include "proto.h"
#include "domain.h"

static FILE *fd;

static void in(int *x, int modus);
static void byten(void *x, int n, int modus);

/* This function reads or writes the restart files.
 * Each processor writes its own restart file, with the
 * I/O being done in parallel. To avoid congestion of the disks
 * you can tell the program to restrict the number of files
 * that are simultaneously written to NumFilesWrittenInParallel.
 *
 * If modus>0  the restart()-routine reads,
 * if modus==0 it writes a restart file.
 */
void restart(int modus)
{
  char buf[200], buf_bak[200], buf_mv[500];
  double save_PartAllocFactor;
  int i, nprocgroup, masterTask, groupTask, old_MaxPart, old_MaxNodes;
  struct global_data_all_processes all_task0;

#if defined(SFR) || defined(BLACK_HOLES)
#ifdef NO_TREEDATA_IN_RESTART
  if(modus == 0)
    {
      rearrange_particle_sequence();
      All.NumForcesSinceLastDomainDecomp = (long long) (1 + All.TreeDomainUpdateFrequency * All.TotNumPart);	/* ensures that new tree will be constructed */
    }
#endif
#endif

  if(ThisTask == 0 && modus == 0)
    {
      sprintf(buf, "%s/restartfiles", All.OutputDir);
      mkdir(buf, 02755);
    }
  MPI_Barrier(MPI_COMM_WORLD);

  sprintf(buf, "%s/restartfiles/%s.%d", All.OutputDir, All.RestartFile, ThisTask);
  sprintf(buf_bak, "%s/restartfiles/%s.%d.bak", All.OutputDir, All.RestartFile, ThisTask);
  sprintf(buf_mv, "mv %s %s", buf, buf_bak);

  if((NTask < All.NumFilesWrittenInParallel))
    {
      printf
	("Fatal error.\nNumber of processors must be a smaller or equal than `NumFilesWrittenInParallel'.\n");
      endrun(2131);
    }

  nprocgroup = NTask / All.NumFilesWrittenInParallel;

  if((NTask % All.NumFilesWrittenInParallel))
    {
      nprocgroup++;
    }

  masterTask = (ThisTask / nprocgroup) * nprocgroup;

  for(groupTask = 0; groupTask < nprocgroup; groupTask++)
    {
      if(ThisTask == (masterTask + groupTask))
	{
	  if(!modus)
	    {
#ifndef NOCALLSOFSYSTEM
	      int ret;

	      ret = system(buf_mv);	/* move old restart files to .bak files */
#endif
	    }
	}
    }

  for(groupTask = 0; groupTask < nprocgroup; groupTask++)
    {
      if(ThisTask == (masterTask + groupTask))	/* ok, it's this processor's turn */
	{
	  if(modus)
	    {
	      if(!(fd = fopen(buf, "r")))
		{
		  printf("Restart file '%s' not found.\n", buf);
		  endrun(7870);
		}
	    }
	  else
	    {
	      if(!(fd = fopen(buf, "w")))
		{
		  printf("Restart file '%s' cannot be opened.\n", buf);
		  endrun(7878);
		}
	    }

	  save_PartAllocFactor = All.PartAllocFactor;

	  /* common data  */
	  byten(&All, sizeof(struct global_data_all_processes), modus);

	  if(ThisTask == 0 && modus > 0)
	    all_task0 = All;

	  if(modus > 0 && groupTask == 0)	/* read */
	    {
	      MPI_Bcast(&all_task0, sizeof(struct global_data_all_processes), MPI_BYTE, 0, MPI_COMM_WORLD);
	    }

	  old_MaxPart = All.MaxPart;
	  old_MaxNodes = (int) (All.TreeAllocFactor * All.MaxPart) + NTopnodes;

	  if(modus)		/* read */
	    {
	      if(All.PartAllocFactor != save_PartAllocFactor)
		{
		  All.PartAllocFactor = save_PartAllocFactor;
		  All.MaxPart = (int) (All.PartAllocFactor * (All.TotNumPart / NTask));
		  All.MaxPartSph = (int) (All.PartAllocFactor * (All.TotN_gas / NTask));
#ifdef INHOMOG_GASDISTR_HINT
		  All.MaxPartSph = All.MaxPart;
#endif
		  save_PartAllocFactor = -1;
		}

	      if(all_task0.Time != All.Time)
		{
		  printf("The restart file on task=%d is not consistent with the one on task=0\n", ThisTask);
		  fflush(stdout);
		  endrun(16);
		}

	      allocate_memory();
	    }

	  in(&NumPart, modus);

	  if(NumPart > All.MaxPart)
	    {
	      printf
		("it seems you have reduced(!) 'PartAllocFactor' below the value of %g needed to load the restart file.\n",
		 NumPart / (((double) All.TotNumPart) / NTask));
	      printf("fatal error\n");
	      endrun(22);
	    }

	  /* Particle data  */
	  byten(&P[0], NumPart * sizeof(struct particle_data), modus);

	  in(&N_gas, modus);

	  if(N_gas > 0)
	    {
	      if(N_gas > All.MaxPartSph)
		{
		  printf
		    ("SPH: it seems you have reduced(!) 'PartAllocFactor' below the value of %g needed to load the restart file.\n",
		     N_gas / (((double) All.TotN_gas) / NTask));
		  printf("fatal error\n");
		  endrun(222);
		}
	      /* Sph-Particle data  */
	      byten(&SphP[0], N_gas * sizeof(struct sph_particle_data), modus);
	    }

	  /* write state of random number generator */
	  byten(gsl_rng_state(random_generator), gsl_rng_size(random_generator), modus);

	  /* now store relevant data for tree */
#ifdef SFR
	  in(&Stars_converted, modus);
#endif

#ifndef NO_TREEDATA_IN_RESTART
	  /* now store relevant data for tree */
	  int nmulti = MULTIPLEDOMAINS;

	  in(&nmulti, modus);
	  if(modus != 0 && nmulti != MULTIPLEDOMAINS)
	    {
	      if(ThisTask == 0)
		printf
		  ("Looks like you changed MULTIPLEDOMAINS from %d to %d.\nWill discard tree stored in restart files and construct a new one.\n",
		   nmulti, (int) MULTIPLEDOMAINS);

	      All.NumForcesSinceLastDomainDecomp = (long long) (1 + All.TreeDomainUpdateFrequency * All.TotNumPart);	/* ensures that new tree will be constructed */
	    }
	  else
	    {
	      in(&NTopleaves, modus);
	      in(&NTopnodes, modus);

	      if(modus)		/* read */
		{
		  domain_allocate();
		  force_treeallocate((int) (All.TreeAllocFactor * All.MaxPart) + NTopnodes, All.MaxPart);
		}

	      in(&Numnodestree, modus);

	      if(Numnodestree > MaxNodes)
		{
		  printf
		    ("Tree storage: it seems you have reduced(!) 'PartAllocFactor' below the value needed to load the restart file (task=%d). "
		     "Numnodestree=%d  MaxNodes=%d\n", ThisTask, Numnodestree, MaxNodes);
		  endrun(221);
		}

	      byten(Nodes_base, Numnodestree * sizeof(struct NODE), modus);
	      byten(Extnodes_base, Numnodestree * sizeof(struct extNODE), modus);

	      byten(Father, NumPart * sizeof(int), modus);

	      byten(Nextnode, NumPart * sizeof(int), modus);
	      byten(Nextnode + All.MaxPart, NTopnodes * sizeof(int), modus);

	      byten(DomainStartList, NTask * MULTIPLEDOMAINS * sizeof(int), modus);
	      byten(DomainEndList, NTask * MULTIPLEDOMAINS * sizeof(int), modus);
	      byten(DomainTask, NTopnodes * sizeof(int), modus);
	      byten(DomainNodeIndex, NTopleaves * sizeof(int), modus);

	      byten(DomainCorner, 3 * sizeof(double), modus);
	      byten(DomainCenter, 3 * sizeof(double), modus);
	      byten(&DomainLen, sizeof(double), modus);
	      byten(&DomainFac, sizeof(double), modus);

	      if(modus)		/* read */
		if(All.PartAllocFactor != save_PartAllocFactor)
		  {
		    for(i = 0; i < NumPart; i++)
		      Father[i] += (All.MaxPart - old_MaxPart);

		    for(i = 0; i < NumPart; i++)
		      if(Nextnode[i] >= old_MaxPart)
			{
			  if(Nextnode[i] >= old_MaxPart + old_MaxNodes)
			    Nextnode[i] += (All.MaxPart - old_MaxPart) + (MaxNodes - old_MaxPart);
			  else
			    Nextnode[i] += (All.MaxPart - old_MaxPart);
			}

		    for(i = 0; i < Numnodestree; i++)
		      {
			if(Nodes_base[i].u.d.sibling >= old_MaxPart)
			  {
			    if(Nodes_base[i].u.d.sibling >= old_MaxPart + old_MaxNodes)
			      Nodes_base[i].u.d.sibling +=
				(All.MaxPart - old_MaxPart) + (MaxNodes - old_MaxNodes);
			    else
			      Nodes_base[i].u.d.sibling += (All.MaxPart - old_MaxPart);
			  }

			if(Nodes_base[i].u.d.father >= old_MaxPart)
			  {
			    if(Nodes_base[i].u.d.father >= old_MaxPart + old_MaxNodes)
			      Nodes_base[i].u.d.father +=
				(All.MaxPart - old_MaxPart) + (MaxNodes - old_MaxNodes);
			    else
			      Nodes_base[i].u.d.father += (All.MaxPart - old_MaxPart);
			  }

			if(Nodes_base[i].u.d.nextnode >= old_MaxPart)
			  {
			    if(Nodes_base[i].u.d.nextnode >= old_MaxPart + old_MaxNodes)
			      Nodes_base[i].u.d.nextnode +=
				(All.MaxPart - old_MaxPart) + (MaxNodes - old_MaxNodes);
			    else
			      Nodes_base[i].u.d.nextnode += (All.MaxPart - old_MaxPart);
			  }
		      }

		    for(i = 0; i < NTopnodes; i++)
		      if(Nextnode[i + All.MaxPart] >= old_MaxPart)
			{
			  if(Nextnode[i + All.MaxPart] >= old_MaxPart + old_MaxNodes)
			    Nextnode[i + All.MaxPart] +=
			      (All.MaxPart - old_MaxPart) + (MaxNodes - old_MaxNodes);
			  else
			    Nextnode[i + All.MaxPart] += (All.MaxPart - old_MaxPart);
			}

		    for(i = 0; i < NTopnodes; i++)
		      if(DomainNodeIndex[i] >= old_MaxPart)
			{
			  if(DomainNodeIndex[i] >= old_MaxPart + old_MaxNodes)
			    DomainNodeIndex[i] += (All.MaxPart - old_MaxPart) + (MaxNodes - old_MaxNodes);
			  else
			    DomainNodeIndex[i] += (All.MaxPart - old_MaxPart);
			}
		  }
	    }
#endif
	  fclose(fd);
	}
      else			/* wait inside the group */
	{
	  if(modus > 0 && groupTask == 0)	/* read */
	    {
	      MPI_Bcast(&all_task0, sizeof(struct global_data_all_processes), MPI_BYTE, 0, MPI_COMM_WORLD);
	    }
	}

      MPI_Barrier(MPI_COMM_WORLD);
    }
#ifdef HEALPIX
  //this should be readed in the parameterfile
  All.Nside=8;
  //
  if(ThisTask==0)
  printf(" Restart calculation of Healpix %i with %i \n",All.Nside,NSIDE2NPIX(All.Nside));
  // initialize the healpix array (just in case)
  All.healpixmap=(float *) malloc(NSIDE2NPIX(All.Nside) * sizeof(float));
  for(i=0;i<NSIDE2NPIX(All.Nside);i++)All.healpixmap[i]=0;
  healpix_halo(All.healpixmap);
#endif
}

/* reads/writes n bytes
 */
void byten(void *x, int n, int modus)
{
  if(modus)
    my_fread(x, n, 1, fd);
  else
    my_fwrite(x, n, 1, fd);
}

/* reads/writes one int
 */
void in(int *x, int modus)
{
  if(modus)
    my_fread(x, 1, sizeof(int), fd);
  else
    my_fwrite(x, 1, sizeof(int), fd);
}

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


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

int main(int argc, char *argv[])
{
    /* 0 + 2x - 3x^2 + 1x^3 */
    double p[] = {0, 2, -3, 1};
    double z[6];
    gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc(4);
    gsl_poly_complex_solve(p, 4, w, z);
    gsl_poly_complex_workspace_free(w);

    for(int i = 0; i < 3; ++i)
        printf("%.12f\n", z[2 * i]);

    return 0;
}

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode