/*<

Transp3d

Usage:
Transp3d  <in.H >out.H 

Input Parameters:

plane       -  int  [12] Plane to transpose
max_memory  -  int  [10] Maximum memory to use in  megabytes
verb        -  int  [0]  Whether or not to be verbose

Output Parameters:



Description:



>*/ 
/*
-------------------------------------------------

Author: Robert Clapp, ESMB 463, 7230253

Date Created:Fri Feb 18 10:11:05 PST 2000

Purpose: 

*/	 

#include <sep.main> 
#include<superset.h>



int read_in[5],rite_in[5],read_out[5],rite_out[5];
char temp_file[1024];
int esize;

#if NeedFunctionPrototypes
_XFUNCPROTOBEGIN
int read_pars(int *verb, int *plane,int *mem, int *nthreads,int *ndim, int *pipe_in, int *pipe_out);
int check_pars(int verb,int *plane,int mem, int nthreads, int ndim);
int rite_pars(int *plane, int mem, int nthreads);
int standardize_it(int verb, int *plane, int mem, int nthreads,int *nin, int *nout,int ndim,int pipein,int pipeout);
int transp_it(int verb,int *nin, int *nout,int mem,int nthreads, int pipein, int pipeout);
int transp_reg(char *tagin,char *tagout,int f5in,int f5out,int nread,int *rite_in,int *rite_out,int *nin,int *nout,char *block_in,char *block_out,int verb);
_XFUNCPROTOEND
#else
int read_pars();
int check_pars();
int rite_pars();
int standardize_it();
int transp_it();
int transp_reg();
#endif



MAIN()
{ 

  int n1;
	int verb,nthreads,mem,ndim;
	int plane[2],pipe_in,pipe_out;
	int nin[5],nout[5];
  char *array;

		seperr("trouble reading parameters \n");
	if(SUCCESS!=check_pars(verb,plane,mem,nthreads,ndim))
		seperr("trouble checking paramters\n");
	if(SUCCESS!=rite_pars(plane,mem,nthreads))
		seperr("trouble writing parameters \n");
	if(SUCCESS!=standardize_it(verb,plane,mem,nthreads,nin,nout,
    ndim,pipe_in,pipe_out))
		seperr("trouble writing parameters \n");

	if(SUCCESS!=transp_it(verb,nin,nout,mem,
   nthreads,pipe_in,pipe_out)){
		seperr("trouble transposing \n");
	}
	
	return(0);
 
 } 
#if NeedFunctionPrototypes
_XFUNCPROTOBEGIN
int read_pars(int *verb,int *plane, int *mem, int *nthreads,int *ndim,int *pipein, int *pipeout)
_XFUNCPROTOEND
#else
int read_pars(verb,plane,mem,nthreads,ndim,pipein,pipeout)
int *verb,*plane,*mem,*nthreads,*ndim,*pipein,*pipeout;
#endif
{
char pl[2],temp_ch[3];

if(0==getch("verb","d",verb)) *verb=0;

*pipein=sseekable("in");
*pipeout=sseekable("out");
if(*verb==1){ 
	fprintf(stderr,"input seekable=%d output seekable=%d \n",*pipein,*pipeout);
}

if(0==getch("nthreads","d",nthreads)) *nthreads=1;

if(0==getch("max_memory","d",mem)) *mem=10;


if(SUCCESS!=sep3d_tag_init("in","in","INPUT"))    
	seperr("trouble initializing tag %s \n","in");

if(SUCCESS!=sep3d_grab_ndims("in",ndim))    
	seperr("trouble  reading ndim from tag in \n");

if(0==getch("plane","s",pl)) strcpy(pl,"12");
	sprintf(temp_ch,"%c",pl[0]); plane[0]=atoi(temp_ch);
	sprintf(temp_ch,"%c",pl[1]); plane[1]=atoi(temp_ch);

;

return(0);


}

#if NeedFunctionPrototypes
_XFUNCPROTOBEGIN
int check_pars(int verb,int *plane, int mem, int nthreads, int ndim)
_XFUNCPROTOEND
#else
int check_pars(verb,plane,mem,nthreads,ndim)
int verb,mem,nthreads,ndim,*plane;
#endif
{
char temp_ch[1024];
int tempi;
float o,d;
char label[128],unit[128];
int n[2];

/* Start with the stupidity checks*/
if(mem< 1) seperr("memory must be 1 or greater \n");

if(nthreads<1) seperr("nthreads must be 1 or greate\n");


/*make first plane smaller axis for convenience*/
if(plane[0] >plane[1]){ tempi=plane[1]; plane[1]=plane[0];plane[0]=tempi;}

if(plane[0] < 1)
	seperr("Bad plane argument, plane=%d data limits 0 to %d \n",plane[0],ndim);
if(plane[1] < 1)
	seperr("Bad plane argument, plane=%d data limits 0 to %d \n",plane[1],ndim);
if(plane[0] == plane[1])
	seperr("Can't transpose a plane with itself \n");

if(SUCCESS!=sep3d_grab_file_type("in",temp_ch))
    seperr("trouble figuring out file type for tag %s \n","in");

if(0!=strcmp(temp_ch,"REGULAR"))
	seperr("Use Sort3d for a non-SEP regular dataset\n");

if(0==hetch("esize","d",&esize)) esize=4;
tempi=esize;

if(SUCCESS!=sep3d_grab_axis("in",plane[0],&n[0],&o,&d,label,unit))
	seperr("trouble reading axis %d from in \n",plane[0]+1);

if(SUCCESS!=sep3d_grab_axis("in",plane[1],&n[1],&o,&d,label,unit))
	seperr("trouble reading axis %d from in \n",plane[1]+1);

if((double)mem*1000000.< 2. * (double) esize * (double)n[0]*(double)n[1])
	seperr("Increase mem, can't hold a single slice in memory\n");

if(verb>0){
	fprintf(stderr,"Transposing plane %d and %d \n",plane[0],plane[1]);
	fprintf(stderr,"Using %d megabytes and %d processors \n",mem,nthreads);
}

return(0);

}

#if NeedFunctionPrototypes
_XFUNCPROTOBEGIN
int rite_pars(int *plane, int mem, int nthreads)
_XFUNCPROTOEND
#else
int rite_pars(plane,mem,nthreads)
int *plane,mem,nthreads;
#endif
{
int i1,ndim;
char label[128],unit[128];
int  n,ierr;
float o,d;
char temp_ch[128];

puthead("Transposing plane %d and %d \n",plane[0],plane[1]);
puthead("Using %d megabytes and %d processors \n",mem,nthreads);


/*create output dataset */


if(SUCCESS!=sep3d_struct_init("in","out","OUTPUT"))
   seperr("trouble creating output structure \n");

if(SUCCESS!=sep3d_grab_axis("in",plane[0],&n,&o,&d,label,unit))
	seperr("trouble reading axis %d from in \n",plane[0]+1);
if(SUCCESS!=sep3d_set_axis("out",plane[1],n,o,d,label,unit))
	seperr("trouble writing axis %d to out \n",plane[1]+1);
if(SUCCESS!=sep3d_grab_axis("in",plane[1],&n,&o,&d,label,unit))
	seperr("trouble reading axis %d from in \n",plane[0]+1);
if(SUCCESS!=sep3d_set_axis("out",plane[0],n,o,d,label,unit))
	seperr("trouble writing axis %d to out \n",plane[0]+1);

if(SUCCESS!=sep3d_rite_format("out","out"))
	seperr("trouble writing out dataset parameters \n");




return(0);
}
#if NeedFunctionPrototypes
_XFUNCPROTOBEGIN
int standardize_it(int verb, int *plane, int mem, int nthreads, int *nin, int *nout,int ndim,int pipein, int pipeout)
_XFUNCPROTOEND
#else
int standardize_it(verb,plane,mem,nthreads,nin,nout,ndim,pipein,pipeout)
int *plane,mem,nthreads,*nin,*nout,ndim,pipein,pipeout;
#endif
{
int n;
float o,d;
double memd;
int axis[6],tempi,i,itemp;
double used,use2,tot[5];
char label[128], unit[128],temp_ch[128];




/*                                                   */
/*  Let us reshape the input and output files        */
/*                                                   */
axis[0]=plane[0]-1;
axis[1]=plane[0];
axis[2]=plane[1]-1;
axis[3]=plane[1];
axis[4]=MAX(ndim,plane[1]);



if(SUCCESS!=sep3d_change_dims("in",5,axis))
	seperr("trouble resetting dims for tag: in \n");
if(SUCCESS!=sep3d_change_dims("out",5,axis))
	seperr("trouble resetting dims for tag: out \n");

for(i=0; i < 5; i++){
	if(SUCCESS!=sep3d_grab_axis("in",i+1,&nin[i],&o,&d,label,unit))
		seperr("trouble grabbing axis %d from tag : in \n",i+1);
	if(SUCCESS!=sep3d_grab_axis("out",i+1,&nout[i],&o,&d,label,unit))
		seperr("trouble grabbing axis %d from tag : out \n",i+1);
}


/*                                                   */
/*   Now lets see if we have to create a temp file   */
/*                                                   */

used=esize;
memd=(double)mem*1000000.;

tot[0]=(double)nin[0]*used; tot[1]=(double)nin[1]*tot[0]; tot[2]=(double)nin[2]*tot[1];
tot[3]=(double)nin[3]*tot[2]; tot[4]=(double)nin[4]*tot[3];


strcpy(temp_file,"in");

if(tot[3] *2. >memd && pipein==1 && pipeout==1){
	if(verb>0) fprintf(stderr,"I am going to have create a temporary file \n");
	itemp=1;
	strcpy(temp_file,"TEMP_XXXXXX"); mktemp(temp_file);
	if(NULL!=auxinout(temp_file))
		seperr("trouble creating temporary file \n");
}
else if(verb>0){
	itemp=0;
	fprintf(stderr,"I don't have to create a temporary file \n");
}

/* The logic is going to be:
     Given a choice loop the output 

*/

if(pipein==1){ /* if the input is a pipe we need to loop the input */
	if(tot[0]>memd) { 
		read_in[0]=MAX(1,(int)(memd/used)); 
		read_in[1]=1; read_in[2]=1; read_in[3]=1; read_in[4]=1;
	}
	else if(tot[1] > memd){
		read_in[0]=nin[0]; read_in[1]=MAX(1,(int)(memd/tot[0])); 
		read_in[2]=1; read_in[3]=1; read_in[4]=1;
	}
	else if(tot[2] > memd){
		read_in[0]=nin[0]; read_in[1]=nin[1]; 
		read_in[2]=MAX(1,(int)(memd/tot[1]));
		read_in[3]=1; read_in[4]=1;
	}
	else if(tot[3]*2 > memd){
		read_in[0]=nin[0]; read_in[1]=nin[1]; 
		read_in[3]=MIN(MAX(1,(int)(memd/tot[2]/2.)),nin[3]);
		read_in[2]=nin[2]; read_in[4]=1;
	}
	else if(tot[4]*2 > memd){
		read_in[0]=nin[0]; read_in[1]=nin[1]; 
		read_in[4]=MAX(1,(int)(memd/tot[3]/2.));
		read_in[2]=nin[2]; read_in[3]=nin[3];
	}
	else{
		read_in[0]=nin[0]; read_in[1]=nin[1]; 
		read_in[4]=nin[4];
		read_in[2]=nin[2]; read_in[3]=nin[3];
	}

	if(pipeout==0){
		read_out[0]=read_in[0]; read_out[1]=read_in[3]; read_out[2]=read_in[2];
		read_out[3]=read_in[1]; read_out[4]=read_in[4];

		rite_in[0]=read_in[0]; rite_in[1]=read_in[1]; rite_in[2]=read_in[2];
		rite_in[3]=read_in[3]; rite_in[4]=read_in[4];

		rite_out[0]=read_in[0]; rite_out[1]=read_in[3]; rite_out[2]=read_in[2];
		rite_out[3]=read_in[1]; rite_out[4]=read_in[4];
	}
	else{
		read_out[0]=read_in[0]; read_out[1]=read_in[1]; read_out[2]=read_in[2];
		read_out[3]=read_in[3]; read_out[4]=read_in[4];
	}
}

tot[0]=(double)nout[0]*used; tot[1]=(double)nout[1]*tot[0]; tot[2]=(double)nout[2]*tot[1];
tot[3]=(double)nout[3]*tot[2]; tot[4]=(double)nout[4]*tot[3];


if(!(pipein==1 && pipeout==0)){ /*if we are looping the output */
	if(tot[0]>memd) { 
		rite_out[0]=MAX(1,(int)(memd/used)); 
		rite_out[1]=1; rite_out[2]=1; rite_out[3]=1; rite_out[4]=1;
	}
	else if(tot[1] > memd){
		rite_out[0]=nout[0]; rite_out[1]=MAX(1,(int)(memd/tot[0])); 
		rite_out[2]=1; rite_out[3]=1; rite_out[4]=1;
	}
	else if(tot[2] > memd){
		rite_out[0]=nout[0]; rite_out[1]=nout[1]; 
		rite_out[2]=MAX(1,(int)(memd/tot[1]));
		rite_out[3]=1; rite_out[4]=1;
	}
	else if(tot[3]*2 > memd){
		rite_out[0]=nout[0]; rite_out[1]=nout[1]; 
		rite_out[3]=MIN(MAX(1,(int)(memd/tot[2]/2.)),nout[3]);
		rite_out[2]=nout[2]; rite_out[4]=1;
	}
	else if(tot[4]*2 > memd){
		rite_out[0]=nout[0]; rite_out[1]=nout[1]; 
		rite_out[4]=MAX(1,(int)(memd/tot[3]/2.));
		rite_out[2]=nout[2]; rite_out[3]=nout[3];
	}
	else{
		rite_out[0]=nout[0]; rite_out[1]=nout[1]; 
		rite_out[4]=nout[4];
		rite_out[2]=nout[2]; rite_out[3]=nout[3];
	}

	rite_in[0]=rite_out[0]; rite_in[1]=rite_out[1]; rite_in[2]=rite_out[2];
	rite_in[3]=rite_out[3]; rite_in[4]=rite_out[4];
	
}
if(verb>0){
	fprintf(stderr,"I am doing my block read_in with %d %d %d %d %d \n",read_in[0],read_in[1],read_in[2],read_in[3],read_in[4]);
	fprintf(stderr,"I am doing my block read_out with %d %d %d %d %d \n",read_out[0],read_out[1],read_out[2],read_out[3],read_out[4]);
	fprintf(stderr,"I am doing my block rite_in with %d %d %d %d %d \n",rite_in[0],rite_in[1],rite_in[2],rite_in[3],rite_in[4]);
	fprintf(stderr,"I am doing my block rite_out with %d %d %d %d %d \n",rite_out[0],rite_out[1],rite_out[2],rite_out[3],rite_out[4]);
}

return(0);
}
#if NeedFunctionPrototypes
_XFUNCPROTOBEGIN
int transp_it(int verb,int *nin, int *nout,int mem,int nthreads, int pipein, int pipeout)
_XFUNCPROTOBEGIN
#else    
int transp_it(verb,nin,nout, mem,nthreads,pipein,pipeout)
int verb,*nin, *nout, mem,nthreads,pipein,pipeout;
#endif
{
int done,nread,fread;
int i4,i3,i2,i1,n,ntot1,ntot2,ntot;
char read_from[128];
char *block_in,*block_out;


	ntot1=rite_out[0]*rite_out[1]*rite_out[2]*rite_out[3]*rite_out[4]*esize;
	if(pipein==1)
		ntot2=read_in[0]*read_in[1]*read_in[2]*read_in[3]*read_in[4]*esize;
	else ntot2=0;

	ntot=MAX(ntot1,ntot2);



	block_in=(char*) alloc(ntot);
	block_out=(char*) alloc(ntot);


done=0;
while(done < nin[4]){
	nread=MIN(MAX(read_in[4],rite_out[4]),nin[4]-done);
	fread=done;

	if(pipein==1 && pipeout==1){
		if(0!=strcmp(temp_file,"in")){
			if(SUCCESS!=transfer_block("in",temp_file,nin,fread,nread,ntot,block_in))
				seperr("trouble copying temp buffer \n");
			if(verb &&  nin[4]>1)
      done+nread, nin[4]);
		}
		if(SUCCESS!=transp_reg(temp_file,"out",0,fread,nread,rite_in,rite_out,
       nin,nout,block_in,block_out,verb))
			seperr("trouble transposing block \n");
	}
	else if(pipeout==0){
		if(SUCCESS!=transp_reg("in","out",fread,fread,nread,rite_in,rite_out,nin,
      nout,block_in,block_out,verb))
			seperr("trouble transposing regular block \n");
	}
	else{
		if(SUCCESS!=transp_reg("in","out",fread,fread,nread,rite_in,rite_out,nin,
      nout,block_in,block_out,verb))
			seperr("trouble transposing regular block \n");
	}
	done=done+nread;
	if(verb && nin[4]>1) fprintf(stderr,"Finished block %d of %d \n",done,nin[4]);
}

return(SUCCESS);
}
#if NeedFunctionPrototypes
_XFUNCPROTOBEGIN
int transp_reg(char *tagin,char *tagout,int f5in,int f5out,int n5,int *rite_in,int *rite_out,int *nin,int *nout,char *block_in,char *block_out,int verb)
_XFUNCPROTOEND
#else    
int transp_reg(tagin,tagout,f5in,f5out,n5,rite_in,rite_out,nin,nout,block_in,block_out,verb)
char *tagin,*tagout;
int f5in, f5out,n5,*rite_in,*rite_out,*nin,*nout,verb;
char *block_in,char *block_out;
#endif
{
int i1,i2,i3,i4,i5;
int dn[4];
int fread[5],frite[5],j[5],nread[5],nrite[5],ipos;
int n4,n3,n2,n1,in5,in4,in3,in2,in1,out1,out2,out3,out,out5,out4,five;
int ibeg,ihold,isend,iend;
float block[100];





fread[4]=f5in; frite[4]=f5out;five=5;
nread[4]=n5; nrite[4]=n5;
j[0]=1; j[1]=1; j[2]=1; j[3]=1; j[4]=1;
frite[0]=0; frite[1]=0; frite[2]=0; frite[3]=0;
while(frite[3] < nout[3]){
	nrite[3]=MIN(rite_out[3],nout[3]-frite[3]);	
	while(frite[2] < nout[2]){
		nrite[2]=MIN(rite_out[2],nout[2]-frite[2]);	
		while(frite[1] < nout[1]){
			nrite[1]=MIN(rite_out[1],nout[1]-frite[1]);	
			while(frite[0] < nout[0]){
				nrite[0]=MIN(rite_out[0],nout[0]-frite[0]);	
				nread[0]=nrite[0];nread[1]=nrite[3];nread[2]=nrite[2];nread[3]=nrite[1];
				fread[0]=frite[0];fread[1]=frite[3];fread[2]=frite[2];fread[3]=frite[1];




				if(0!=sreed_window(tagin,&five,nin,nread,fread,j,esize,block_in))
					seperr("trouble reading block \n");

					n4=nrite[3];n3=nrite[2];n2=nrite[1];n1=nrite[0];in1=esize;out1=esize;
					in2=n1*in1; out2=n1*out1; 
					in3=in2*n4; out3=out2*n2; 
					in4=in3*n3;out4=out3*n3;
					in5=in4*n2;out5=out4*n4;
					ipos=0;

					
					ibeg=0;
					/*can you believe all of this crap just do do this !!!!*/
					for(i5=ibeg; i5 < n5; i5++){
					for(i4=0; i4 < n4; i4++){
						for(i3=0; i3 < n3; i3++){
							for(i2=0; i2 < n2; i2++){
									memcpy( (void*)(block_out+i2*out2+i3*out3+i4*out4+out5*i5),
										(const void*)(block_in+i4*in2+i3*in3+i2*in4+in5*i5),
                 in2);
								  
							}
						}
					}}

					if(0!=srite_window(tagout,&five,nout,nrite,frite,j,esize,block_out))
						seperr("trouble reading block \n");
				frite[0]=frite[0]+nrite[0]; 
			}
			frite[1]=frite[1]+nrite[1];  frite[0]=0;
		}
		frite[2]=frite[2]+nrite[2]; frite[1]=0;
	}
	frite[3]=frite[3]+nrite[3]; frite[2]=0;
	if(nout[4]==1 && verb) fprintf(stderr,"Finished block %d of %d \n",frite[3],nout[3]);
}

return(SUCCESS);
}
#if NeedFunctionPrototypes
_XFUNCPROTOBEGIN
int transfer_block(char *tagin,char *tagout, int *nin, int f5, int n5, int ntot, char *block_in)
_XFUNCPROTOEND
#else    
int transfer_block(tagin,tagout, nin, f5, n5,ntot,block_in)
char *tagin, *tagout,*block_in; 
int *nin, f5,  n5,ntot;
#endif    
{
double nbig;
double ndone;
int block;



if(f5!=0) sseek(tagout,0,0);
nbig=(double)nin[0]*(double)nin[1]*(double)nin[2]*(double)nin[3];
ndone=0;

while(ndone< nbig){
	block=MIN((int)(nbig-ndone),ntot);
	if(block*esize!=sreed(tagin,block_in,esize*block))
		seperr("troulbe reading from tag=%s \n",tagin);
	if(block*esize!=srite(tagout,block_in,esize*block))
		seperr("troulbe writing to tag=%s \n",tagout);

	ndone=ndone+block;
}


return(SUCCESS);
}