[Genabel-commits] r1469 - in pkg: . OmicABELnoMM
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Dec 16 09:37:16 CET 2013
Author: lckarssen
Date: 2013-12-16 09:37:15 +0100 (Mon, 16 Dec 2013)
New Revision: 1469
Added:
pkg/OmicABELnoMM/
pkg/OmicABELnoMM/AIOwrapper.cpp
pkg/OmicABELnoMM/AIOwrapper.h
pkg/OmicABELnoMM/Algorithm.cpp
pkg/OmicABELnoMM/Algorithm.h
pkg/OmicABELnoMM/CreateData.R
pkg/OmicABELnoMM/Definitions.h
pkg/OmicABELnoMM/Utility.cpp
pkg/OmicABELnoMM/Utility.h
pkg/OmicABELnoMM/XL.fvd
pkg/OmicABELnoMM/XL.fvi
pkg/OmicABELnoMM/XR.fvd
pkg/OmicABELnoMM/XR.fvi
pkg/OmicABELnoMM/Y.fvd
pkg/OmicABELnoMM/Y.fvi
pkg/OmicABELnoMM/howtocompile.txt
pkg/OmicABELnoMM/main.cpp
Log:
First commit of OmicABELnoMM, written by Alvaro Jesus Frank. The name refers to OmicABEL, but without support for mixed models.
Note that OmicABELnoMM is under full development.
Added: pkg/OmicABELnoMM/AIOwrapper.cpp
===================================================================
--- pkg/OmicABELnoMM/AIOwrapper.cpp (rev 0)
+++ pkg/OmicABELnoMM/AIOwrapper.cpp 2013-12-16 08:37:15 UTC (rev 1469)
@@ -0,0 +1,778 @@
+#include "AIOwrapper.h"
+
+
+
+AIOwrapper::AIOwrapper()
+{
+ Fhandler = &FHandler;
+ io_overhead = "*";
+}
+
+AIOwrapper::~AIOwrapper()
+{
+
+}
+
+
+void AIOwrapper::initialize(struct Settings ¶ms)
+{
+
+
+ pthread_mutex_init( &(FHandler.m_more), NULL);
+ pthread_mutex_init( &(FHandler.m_read), NULL);
+ pthread_mutex_init( &(FHandler.m_buff_upd), NULL);
+
+ pthread_cond_init(&(FHandler.condition_more), NULL);
+ pthread_cond_init(&(FHandler.condition_read), NULL);
+
+ pthread_barrier_init(&(FHandler.finalize_barrier),NULL,2);
+
+
+ Fhandler->fakefiles = params.use_fake_files;
+
+ databel_fvi* Yfvi;
+ databel_fvi* ALfvi;
+ databel_fvi* ARfvi;
+
+ if(!Fhandler->fakefiles)
+ {
+ Fhandler->fnameAL = params.fnameAL;
+ Fhandler->fnameAR = params.fnameAR;
+ Fhandler->fnameY = params.fnameY;
+ Fhandler->fnameOutB = params.fnameOutB;
+
+ databel_fvi* Yfvi = load_databel_fvi( (Fhandler->fnameY+".fvi").c_str() );
+ databel_fvi* ALfvi = load_databel_fvi( (Fhandler->fnameAL+".fvi").c_str() );
+ databel_fvi* ARfvi = load_databel_fvi( (Fhandler->fnameAR+".fvi").c_str() );
+ params.n = ALfvi->fvi_header.numObservations;
+ params.m = ARfvi->fvi_header.numVariables/params.r;
+ params.t = Yfvi->fvi_header.numVariables;
+ params.l = ALfvi->fvi_header.numVariables;
+
+ //block size to keep mem under 1 gigabyte
+ int opt_block = params.n/(4*1000^3)*(1/(2*params.r));
+ int opt_tb = 1000;
+ int opt_mb = 1000;
+
+ params.mb = min(params.m,opt_tb);
+ params.tb = min(params.t,opt_mb);
+ }
+ else
+ {
+
+ }
+
+ params.p = params.l + params.r;
+
+
+ //block size to keep mem under 1 gigabyte
+// int opt_block = params.n/(4*1000^3)*(1/(2*params.r));
+// int opt_tb = max(4*2000,opt_block);
+// int opt_mb = max(2000,opt_block);
+//
+// params.mb = min(params.m,opt_tb);
+// params.tb = min(params.t,opt_mb);
+
+ prepare_B();
+ prepare_AL(params.l,params.n);
+ prepare_AR( params.mb, params.n, params.m, params.r);
+ prepare_Y(params.tb, params.n, params.t);
+
+
+
+
+
+}
+
+void AIOwrapper::finalize()
+{
+ //cout << "f";
+ void *status;
+
+ Fhandler->not_done = false;
+ pthread_mutex_lock(&(Fhandler->m_more));
+ pthread_cond_signal( &(Fhandler->condition_more ));
+ pthread_mutex_unlock(&(Fhandler->m_more));
+
+ pthread_barrier_wait(&(Fhandler->finalize_barrier));
+ finalize_Y();
+ finalize_AR();
+ finalize_AL();
+ finalize_B();
+
+ pthread_attr_destroy(&(Fhandler->attr));
+
+ pthread_mutex_destroy(&(Fhandler->m_more));
+ pthread_cond_destroy(&(Fhandler->condition_more));
+
+ pthread_mutex_destroy(&(Fhandler->m_read));
+ pthread_cond_destroy(&(Fhandler->condition_read));
+
+
+}
+
+
+void AIOwrapper::prepare_B()
+{
+ Fhandler->fp_B = fopen((Fhandler->fnameOutB+".fvd").c_str(), "w+b");
+ if(Fhandler->fp_B == 0)
+ {
+ cout << "Error Opening File B " << Fhandler->fnameOutB << endl;
+ exit(1);
+ }
+}
+
+
+void AIOwrapper::finalize_B()
+{
+ fclose(Fhandler->fp_B);
+}
+
+
+void* AIOwrapper::async_io( void *ptr )
+{
+ //cout << "async_io\n" << flush;
+ type_fileh* Fhandler = (type_fileh *)ptr;
+ int size_buff,tmp_y_blockSize,tmp_ar_blockSize;
+
+ struct timespec timeToWait;
+ FILE* fp_Y;
+
+ FILE* fp_Ar;
+ if(!Fhandler->fakefiles)
+ {
+ fp_Y = fopen((Fhandler->fnameY+".fvd").c_str(), "rb");
+ if(fp_Y == 0)
+ {
+ cout << "Error Reading File Y " << Fhandler->fnameY << endl;
+ exit(1);
+ }
+
+ fp_Ar = fopen((Fhandler->fnameAR+".fvd").c_str(), "rb");
+ if(fp_Ar == 0)
+ {
+ cout << "Error Reading File Xr " << Fhandler->fnameAR << endl;
+ exit(1);
+ }
+
+
+ }
+ else
+ {
+// fp_Y = fopen("tempY.bin", "w+b");
+// if(fp_Y == 0)
+// {
+// cout << "Error creating temp File Y " << Fhandler->fnameY << endl;
+// exit(1);
+// }
+// fwrite(Fhandler->Yb, sizeof(type_precision), Fhandler->n*Fhandler->y_blockSize, fp_Y);
+// fclose(fp_Y);
+//
+//
+// fp_Ar = fopen("tempAR.bin", "w+b");
+// if(fp_Ar == 0)
+// {
+// cout << "Error creating temp File Y " << Fhandler->fnameY << endl;
+// exit(1);
+// }
+// fwrite(Fhandler->Yb, sizeof(type_precision), Fhandler->n*Fhandler->y_blockSize, fp_Ar);
+// fclose(fp_Ar);
+//
+// fp_Y = fopen("tempY.bin", "rb");
+// if(fp_Y == 0)
+// {
+// cout << "Error Reading File Y " << Fhandler->fnameY << endl;
+// exit(1);
+// }
+//
+// fp_Ar = fopen("tempAR.bin", "rb");
+// if(fp_Ar == 0)
+// {
+// cout << "Error Reading File Xr " << Fhandler->fnameAR << endl;
+// exit(1);
+// }
+
+ }
+
+
+ Fhandler->not_done = true;
+ Fhandler->reset_wait = false;
+
+ //cout << "c" << flush;
+
+ while(Fhandler->not_done)
+ {
+
+ while(!Fhandler->empty_buffers.empty() && Fhandler->y_to_readSize)
+ {
+
+ tmp_y_blockSize = Fhandler->y_blockSize;
+ if(Fhandler->y_to_readSize < Fhandler->y_blockSize)
+ tmp_y_blockSize = Fhandler->y_to_readSize;
+
+ Fhandler->y_to_readSize -= tmp_y_blockSize;
+ size_buff = Fhandler->n * tmp_y_blockSize;
+ //cout << Fhandler->y_to_readSize << endl;
+
+ pthread_mutex_lock(&(Fhandler->m_buff_upd));
+ //cout << " pre;" << Fhandler->full_buffers.size() << ";" << Fhandler->empty_buffers.size() << endl;
+ type_buffElement* tobeFilled = Fhandler->empty_buffers.front();
+ Fhandler->empty_buffers.pop();
+ //pthread_mutex_unlock(&(Fhandler->m_buff_upd));
+
+ tobeFilled->size = tmp_y_blockSize;
+
+ if(Fhandler->fakefiles)
+ {
+ //fread (tobeFilled->buff,sizeof(type_precision),size,fp_Y);
+ //fseek ( fp_Y , 0 , SEEK_SET );
+ int old_seed = Fhandler->seed;
+ srand (old_seed);
+ re_random_vec(tobeFilled->buff, size_buff );
+ re_random_vec_nan(tobeFilled->buff, size_buff );
+ Fhandler->seed += 75;
+ }
+ else
+ {
+ fread (tobeFilled->buff,sizeof(type_precision),size_buff,fp_Y);
+ if(Fhandler->y_to_readSize <= 0)
+ {
+ fseek ( fp_Y , 0 , SEEK_SET );
+ }
+ }
+
+
+ //pthread_mutex_lock(&(Fhandler->m_buff_upd));
+ Fhandler->full_buffers.push(tobeFilled);
+ // cout << "\nStoring " << tobeFilled << endl;
+ //cout << " post;" << Fhandler->full_buffers.size() << ";" << Fhandler->empty_buffers.size() << endl;
+ pthread_mutex_unlock(&(Fhandler->m_buff_upd));
+
+ pthread_mutex_lock(&(Fhandler->m_read));
+ pthread_cond_signal( &(Fhandler->condition_read ));
+ pthread_mutex_unlock(&(Fhandler->m_read));
+
+ }
+
+ while(!Fhandler->ar_empty_buffers.empty() && Fhandler->Ar_to_readSize)
+ {
+ tmp_ar_blockSize = Fhandler->Ar_blockSize;
+ if(Fhandler->Ar_to_readSize < Fhandler->Ar_blockSize)
+ tmp_ar_blockSize = Fhandler->Ar_to_readSize;
+
+ Fhandler->Ar_to_readSize -= tmp_ar_blockSize;
+ size_buff = Fhandler->n * tmp_ar_blockSize*Fhandler->r;
+
+ pthread_mutex_lock(&(Fhandler->m_buff_upd));
+ type_buffElement* tobeFilled = Fhandler->ar_empty_buffers.front();
+ Fhandler->ar_empty_buffers.pop();
+
+
+ tobeFilled->size = tmp_ar_blockSize;
+
+ if(Fhandler->fakefiles)
+ {
+ //fread (tobeFilled->buff,sizeof(type_precision),size,fp_Ar);
+ //fseek ( fp_Ar , 0 , SEEK_SET );
+ re_random_vec(tobeFilled->buff , Fhandler->n * tmp_ar_blockSize*Fhandler->r );
+ re_random_vec_nan(tobeFilled->buff , Fhandler->n * tmp_ar_blockSize*Fhandler->r );
+
+ }
+ else
+ {
+ fread (tobeFilled->buff,sizeof(type_precision),size_buff,fp_Ar);
+ if(Fhandler->Ar_to_readSize <= 0)
+ {
+ fseek ( fp_Ar , 0 , SEEK_SET );
+ }
+ }
+
+ Fhandler->ar_full_buffers.push(tobeFilled);
+ // cout << "\nStoring " << tobeFilled << endl;
+ pthread_mutex_unlock(&(Fhandler->m_buff_upd));
+
+ pthread_mutex_lock(&(Fhandler->m_read));
+ pthread_cond_signal( &(Fhandler->condition_read ));
+ pthread_mutex_unlock(&(Fhandler->m_read));
+
+ }
+
+
+
+
+#ifdef WINDOWS
+ SYSTEMTIME time;
+ GetSystemTime(&time);
+
+ timeToWait.tv_sec = time.wSecond + 500/1000;
+ long int morenanos = (500%1000)*1000000;
+ timeToWait.tv_nsec = time.wMilliseconds*1000 + morenanos ;
+#else
+ clock_gettime(CLOCK_REALTIME, &timeToWait);
+ timeToWait.tv_nsec += 150;
+#endif
+
+ pthread_mutex_lock(&(Fhandler->m_more));
+ pthread_cond_timedwait( &(Fhandler->condition_more), &(Fhandler->m_more), &timeToWait );
+ pthread_mutex_unlock( &(Fhandler->m_more ));
+
+ pthread_mutex_lock(&(Fhandler->m_read));
+ pthread_cond_signal( &(Fhandler->condition_read ));
+ pthread_mutex_unlock(&(Fhandler->m_read));
+
+ if(Fhandler->reset_wait)
+ {
+ pthread_barrier_wait(&(Fhandler->finalize_barrier));
+ //wait for main thread to reset everything
+ pthread_barrier_wait(&(Fhandler->finalize_barrier));
+ }
+
+
+ }
+ //cout << "k" << flush;
+ //barrier
+ pthread_barrier_wait(&(Fhandler->finalize_barrier));
+
+ type_buffElement* tmp;
+ while(!Fhandler->full_buffers.empty())
+ {
+ tmp= Fhandler->full_buffers.front();
+ Fhandler->full_buffers.pop();
+ delete []tmp->buff;
+ delete tmp;
+ }
+
+ while(!Fhandler->empty_buffers.empty())
+ {
+ tmp= Fhandler->empty_buffers.front();
+ Fhandler->empty_buffers.pop();
+ delete []tmp->buff;
+ delete tmp;
+ }
+
+ while(!Fhandler->ar_full_buffers.empty())
+ {
+ tmp= Fhandler->ar_full_buffers.front();
+ Fhandler->ar_full_buffers.pop();
+ delete tmp->buff;
+ delete tmp;
+ }
+
+ while(!Fhandler->ar_empty_buffers.empty())
+ {
+ tmp= Fhandler->ar_empty_buffers.front();
+ Fhandler->ar_empty_buffers.pop();
+ delete []tmp->buff;
+ delete tmp;
+ }
+
+ //
+
+ //pthread_exit(NULL);
+
+ if(!Fhandler->fakefiles)
+ {
+ fclose(fp_Y);
+ fclose(fp_Ar);
+
+ }
+ else
+ {
+ fclose(fp_Y);
+ fclose(fp_Ar);
+ }
+
+//
+// //!induce realistic fileread delay
+
+}
+
+void AIOwrapper::load_ARblock(type_precision** Ar, int &Ar_blockSize)
+{
+
+ int status;
+ int createstatus = 0;
+ //cout<<"^";
+
+ while(Fhandler->ar_full_buffers.empty())
+ {
+ pthread_mutex_lock(&(Fhandler->m_more));
+ pthread_cond_signal( &(Fhandler->condition_more ));
+ pthread_mutex_unlock(&(Fhandler->m_more));
+
+ io_overhead = "#";
+
+ pthread_mutex_lock(&(Fhandler->m_read));
+ pthread_cond_wait( &(Fhandler->condition_read), &(Fhandler->m_read ));
+ pthread_mutex_unlock(&(Fhandler->m_read));
+ }
+
+
+ //!read new rdy buffer
+ pthread_mutex_lock(&(Fhandler->m_buff_upd));
+ if(Fhandler->Ar_currentReadBuff)
+ {
+ Fhandler->ar_empty_buffers.push(Fhandler->Ar_currentReadBuff);
+ }
+
+ Fhandler->Ar_currentReadBuff = Fhandler->ar_full_buffers.front();
+ Fhandler->ar_full_buffers.pop();
+
+ //cout << "\nReading " << Fhandler->Ar_currentReadBuff << endl;
+ Fhandler->Ar = Fhandler->Ar_currentReadBuff->buff;
+ Ar_blockSize = Fhandler->Ar_currentReadBuff->size;
+
+ pthread_mutex_unlock(&(Fhandler->m_buff_upd));
+
+
+
+ pthread_mutex_lock(&(Fhandler->m_more));
+ pthread_cond_signal( &(Fhandler->condition_more ));
+ pthread_mutex_unlock(&(Fhandler->m_more));
+
+
+
+//secueantial
+// int tmp_ar_blockSize = Fhandler->Ar_blockSize;
+// if(Fhandler->Ar_to_readSize < Fhandler->Ar_blockSize)
+// tmp_ar_blockSize = Fhandler->Ar_to_readSize;
+//
+// Fhandler->Ar_to_readSize -= tmp_ar_blockSize;
+// int size_buff = Fhandler->n * tmp_ar_blockSize * Fhandler->r;
+// re_random_vec(Fhandler->Ar,size_buff);
+//Ar_blockSize=tmp_ar_blockSize;
+
+
+ (*Ar) = Fhandler->Ar;
+
+
+
+
+
+
+}
+
+void AIOwrapper::load_Yblock(type_precision** Y, int &y_blockSize)
+{
+
+ int status;
+ int createstatus = 0;
+
+ while(Fhandler->full_buffers.empty())
+ {
+
+ pthread_mutex_lock(&(Fhandler->m_more));
+ pthread_cond_signal( &(Fhandler->condition_more ));
+ pthread_mutex_unlock(&(Fhandler->m_more));
+
+ io_overhead = "!";
+
+ pthread_mutex_lock(&(Fhandler->m_read));
+ pthread_cond_wait( &(Fhandler->condition_read), &(Fhandler->m_read ));
+ pthread_mutex_unlock(&(Fhandler->m_read));
+
+ }
+
+ //!read new rdy buffer
+ pthread_mutex_lock(&(Fhandler->m_buff_upd));
+ //cout << " pre," << Fhandler->full_buffers.size() << ";" << Fhandler->empty_buffers.size() << endl;
+
+ if(Fhandler->currentReadBuff)
+ {
+ //memset(Fhandler->currentReadBuff->buff,0,y_blockSize);
+ Fhandler->empty_buffers.push(Fhandler->currentReadBuff);
+ }
+ Fhandler->currentReadBuff = Fhandler->full_buffers.front();
+ Fhandler->full_buffers.pop();
+
+ //cout << "\nReading " << Fhandler->currentReadBuff << endl;
+ Fhandler->Yb = Fhandler->currentReadBuff->buff;
+ y_blockSize = Fhandler->currentReadBuff->size;
+
+ (*Y) = Fhandler->Yb;
+
+ //cout << " post," << Fhandler->full_buffers.size() << ";" << Fhandler->empty_buffers.size() << endl;
+
+ pthread_mutex_unlock(&(Fhandler->m_buff_upd));
+
+
+
+ pthread_mutex_lock(&(Fhandler->m_more));
+ pthread_cond_signal( &(Fhandler->condition_more ));
+ pthread_mutex_unlock(&(Fhandler->m_more));
+
+ //matlab_print_matrix("Y",Fhandler->n,y_blockSize,*Y);
+
+
+}
+
+void AIOwrapper::write_B(type_precision* B, int p, int blockSize)
+{
+ fwrite (B , sizeof(type_precision), p*blockSize, Fhandler->fp_B);
+}
+
+void AIOwrapper::prepare_Y(int y_blockSize, int n, int totalY)
+{
+ //for fake files
+
+ Fhandler->seed = 1337;
+ srand (Fhandler->seed);
+
+ Fhandler->y_blockSize = y_blockSize;
+
+ Fhandler->n= n;
+ Fhandler->Y_Amount=totalY;
+ Fhandler->y_to_readSize = Fhandler->Y_Amount;
+ Fhandler->buff_count = min(2,(totalY+ y_blockSize - 1)/y_blockSize) ;
+ //cout << "buffcount " << Fhandler->buff_count;
+
+
+ Fhandler->currentReadBuff = 0;
+ type_buffElement* tmp;
+
+ for(int i = 0; i< Fhandler->buff_count ; i++)
+ {
+ tmp = new type_buffElement();
+ tmp->buff = new type_precision[Fhandler->n*Fhandler->y_blockSize];
+ tmp->size = y_blockSize;
+// for( int i = 0; i < Fhandler->n*Fhandler->y_blockSize; i++)
+// {
+// (tmp->buff)[i] = 0;
+// }
+ Fhandler->empty_buffers.push(tmp);
+ }
+ Fhandler->Yb = tmp->buff;
+
+
+
+ pthread_mutex_init(&(Fhandler->m_buff_upd), NULL);
+ pthread_mutex_init(&(Fhandler->m_more), NULL);
+ pthread_mutex_init(&(Fhandler->m_read), NULL);
+
+ pthread_attr_init(&(Fhandler->attr));
+ pthread_attr_setdetachstate(&(Fhandler->attr), PTHREAD_CREATE_JOINABLE);
+
+ pthread_create( &(Fhandler->iothread),&(Fhandler->attr), AIOwrapper::async_io, (void*)Fhandler);
+
+
+
+}
+
+void AIOwrapper::reset_Y()
+{
+ void *status;
+
+ Fhandler->seed = 1337;
+
+ //cout << "ry" << flush;
+
+ Fhandler->reset_wait = true;
+ pthread_barrier_wait(&(Fhandler->finalize_barrier));
+
+ pthread_mutex_lock(&(Fhandler->m_buff_upd));
+ Fhandler->y_to_readSize = Fhandler->Y_Amount;
+
+ if(Fhandler->currentReadBuff)
+ {
+ Fhandler->full_buffers.push(Fhandler->currentReadBuff);
+ Fhandler->currentReadBuff=0;
+ }
+
+ while(!Fhandler->full_buffers.empty())
+ {
+ Fhandler->empty_buffers.push(Fhandler->full_buffers.front());
+ for( int i = 0; i < Fhandler->n*Fhandler->y_blockSize; i++)
+ {
+ ((Fhandler->full_buffers.front())->buff)[i] = 0;
+ }
+ Fhandler->full_buffers.pop();
+ }
+ pthread_mutex_unlock(&(Fhandler->m_buff_upd));
+
+ Fhandler->reset_wait = false;
+
+ pthread_barrier_wait(&(Fhandler->finalize_barrier));
+
+ pthread_mutex_lock(&(Fhandler->m_more));
+ pthread_cond_signal( &(Fhandler->condition_more ));
+ pthread_mutex_unlock(&(Fhandler->m_more));
+
+
+}
+
+void AIOwrapper::reset_AR()
+{
+ void *status;
+
+
+ //cout << "ra" << flush;
+
+ Fhandler->reset_wait = true;
+ pthread_barrier_wait(&(Fhandler->finalize_barrier));
+
+ pthread_mutex_lock(&(Fhandler->m_buff_upd));
+ Fhandler->Ar_to_readSize = Fhandler->Ar_Amount;
+
+ if(Fhandler->Ar_currentReadBuff)
+ {
+ Fhandler->ar_full_buffers.push(Fhandler->Ar_currentReadBuff);
+ Fhandler->Ar_currentReadBuff=0;
+ }
+
+ while(!Fhandler->ar_full_buffers.empty())
+ {
+ Fhandler->ar_empty_buffers.push(Fhandler->ar_full_buffers.front());
+// for( int i = 0; i < Fhandler->n*Fhandler->r*Fhandler->Ar_blockSize; i++)
+// {
+// ((Fhandler->ar_full_buffers.front())->buff)[i] = 0;
+// }
+ Fhandler->ar_full_buffers.pop();
+ }
+ pthread_mutex_unlock(&(Fhandler->m_buff_upd));
+
+ Fhandler->reset_wait = false;
+
+ pthread_barrier_wait(&(Fhandler->finalize_barrier));
+
+ pthread_mutex_lock(&(Fhandler->m_more));
+ pthread_cond_signal( &(Fhandler->condition_more ));
+ pthread_mutex_unlock(&(Fhandler->m_more));
+
+
+}
+
+void AIOwrapper::finalize_Y()
+{
+
+
+}
+
+void AIOwrapper::prepare_AR( int desired_blockSize, int n, int totalR, int columnsAR)
+{
+
+ Fhandler->Ar = new type_precision[desired_blockSize*columnsAR*n];
+ Fhandler->Ar_blockSize = desired_blockSize;
+ Fhandler->r = columnsAR;
+ Fhandler->Ar_Amount = totalR;
+ Fhandler->Ar_to_readSize = Fhandler->Ar_Amount;
+
+ int buff_count = min(2,(totalR+ desired_blockSize - 1)/desired_blockSize);
+
+ Fhandler->Ar_currentReadBuff = 0;
+ type_buffElement* tmp;
+
+ for(int i = 0; i< buff_count ; i++)
+ {
+ tmp = new type_buffElement();
+ tmp->buff = new type_precision[n*desired_blockSize*columnsAR];
+ tmp->size = desired_blockSize;
+// memset(tmp->buff,-9,n*desired_blockSize*columnsAR);
+
+ Fhandler->ar_empty_buffers.push(tmp);
+ }
+ Fhandler->Ar = tmp->buff;
+
+}
+
+
+void AIOwrapper::finalize_AR()
+{
+
+}
+
+
+void AIOwrapper::load_AL(type_precision** AL)
+{
+ if(Fhandler->fakefiles)
+ {
+ re_random_vec(Fhandler->AL,Fhandler->n*Fhandler->l);
+ re_random_vec_nan(Fhandler->AL,Fhandler->n*Fhandler->l);
+ (*AL) = Fhandler->AL;
+ }
+ else
+ {
+ FILE *fp;
+ fp = fopen((Fhandler->fnameAL+".fvd").c_str(), "rb");
+ if(fp == 0)
+ {
+ cout << "Error Reading File " << Fhandler->fnameAL << endl;
+ exit(1);
+ }
+
+ fread (Fhandler->AL,sizeof(type_precision),Fhandler->l*Fhandler->n,fp);
+ fclose(fp);
+ }
+
+
+
+ (*AL) = Fhandler->AL;
+
+ //matlab_print_matrix("AL",Fhandler->n,Fhandler->l,*AL);
+}
+
+void AIOwrapper::prepare_AL( int columnsAL, int n)
+{
+ Fhandler->AL = new type_precision[columnsAL*n];
+ Fhandler->l=columnsAL;
+}
+
+void AIOwrapper::finalize_AL()
+{
+ delete []Fhandler->AL;
+}
+
+void AIOwrapper::free_databel_fvi( struct databel_fvi **fvi )
+{
+ free ((*fvi)->fvi_data);
+ free (*fvi);
+ *fvi = NULL;
+}
+
+FILE * AIOwrapper::fgls_fopen( const char *path, const char *mode )
+{
+ FILE * f;
+ char err[100];
+
+ f = fopen( path, mode );
+ if ( f == NULL )
+ {
+ cout << "\nerror on fgls_fopen\n";
+ exit( 1 );
+ }
+ return f;
+}
+
+void * AIOwrapper::fgls_malloc_impl( const char* file, long line, size_t size )
+{
+ void *buf;
+
+ if ( (buf = malloc( size )) == NULL ) {
+ cout<< "\nCouldn't allocate %ld bytes of memory in %s:%ld\n";
+ exit(1);
+ }
+
+ return buf;
+}
+
+struct databel_fvi * AIOwrapper::load_databel_fvi( const char *path )
+{
+ FILE *f;
+ databel_fvi *fvi;
+ size_t data_size;
+
+ f = fgls_fopen( path, "r" );
+
+ fvi = (databel_fvi*) fgls_malloc( sizeof(databel_fvi) );
+ // Header
+ fread( &fvi->fvi_header, sizeof(databel_fvi_header), 1, f );
+ // Labels
+ data_size = (fvi->fvi_header.numVariables +fvi->fvi_header.numObservations ) *
+ fvi->fvi_header.namelength * sizeof(char);
+ fvi->fvi_data = (char *) fgls_malloc ( data_size );
+ // Load labels
+ fread( fvi->fvi_data, 1, data_size, f );
+
+ fclose( f );
+
+ return fvi;
+}
Added: pkg/OmicABELnoMM/AIOwrapper.h
===================================================================
--- pkg/OmicABELnoMM/AIOwrapper.h (rev 0)
+++ pkg/OmicABELnoMM/AIOwrapper.h 2013-12-16 08:37:15 UTC (rev 1469)
@@ -0,0 +1,178 @@
+#ifndef AIOWRAPPER_H
+#define AIOWRAPPER_H
+
+#include "Definitions.h"
+#include "Utility.h"
+
+typedef struct BufferElement type_buffElement;
+
+struct BufferElement
+{
+ type_precision* buff;
+ int size;
+};
+
+typedef struct fileh type_fileh;
+
+
+
+struct fileh
+{
+ string fnameAL;
+ string fnameAR;
+ string fnameY;
+ string fnameOutB;
+
+ FILE* fp_B;
+
+ bool doublefileType;
+ bool fakefiles;
+
+ type_precision* Yb;
+ type_precision* Ar;
+ type_precision* AL;
+ type_buffElement* currentReadBuff;
+ type_buffElement* Ar_currentReadBuff;
+ int buff_count;
+
+ queue<type_buffElement*> empty_buffers;
+ queue<type_buffElement*> full_buffers;
+
+ queue<type_buffElement*> b_empty_buffers;
+ queue<type_buffElement*> b_full_buffers;
+
+ queue<type_buffElement*> ar_empty_buffers;
+ queue<type_buffElement*> ar_full_buffers;
+
+ int index;
+
+ int n;
+ int r;
+ int l;
+
+ int Ar_Amount;
+ int Ar_blockSize;
+ int Ar_to_readSize;
+
+ int Y_Amount;
+ int y_blockSize;
+ int y_to_readSize;
+
+ bool not_done;
+ bool reset_wait;
+
+ int seed;
+ int Aseed;
+
+ pthread_mutex_t m_more ;
+ pthread_cond_t condition_more ;
+ pthread_mutex_t m_read ;
+
+ pthread_cond_t condition_read ;
+ pthread_mutex_t m_buff_upd ;
+
+ pthread_t iothread;
+ pthread_attr_t attr;
+
+ //barrier
+ pthread_barrier_t finalize_barrier;
+
+
+};
+
+#define fgls_malloc(size) fgls_malloc_impl(__FILE__, __LINE__, size)
+
+enum datatype{ UNSIGNED_SHORT_INT_TYPE = 1,
+ SHORT_INT_TYPE,
+ UNSIGNED_INT_TYPE,
+ INT_TYPE,
+ FLOAT_TYPE,
+ DOUBLE_TYPE,
+ SIGNED_CHAR_TYPE,
+ UNSIGNED_CHAR_TYPE };
+
+#define NAMELENGTH 32
+#define RESERVEDSPACE 5
+
+typedef struct databel_fvi_header
+{
+ unsigned short int type;
+ unsigned int nelements;
+ unsigned int numObservations;
+ unsigned int numVariables;
+ unsigned int bytesPerRecord;
+ unsigned int bitsPerRecord;
+ unsigned int namelength;
+ unsigned int reserved[RESERVEDSPACE];
+} databel_fvi_header;
+
+typedef struct databel_fvi
+{
+ databel_fvi_header fvi_header;
+ char *fvi_data;
+} databel_fvi;
+
+class AIOwrapper
+{
+ public:
+ AIOwrapper();
+ ~AIOwrapper();
+
+ void initialize(struct Settings ¶ms);
+ void finalize();
+
+
+ void load_AL(type_precision** AL);
+ void load_ARblock(type_precision** Y, int &blockSize);
+ void load_Yblock(type_precision** Y, int &blockSize);
+ void reset_Y();
+ void reset_AR();
+
+ void write_B(type_precision* B, int p, int blockSize);
+
+ string io_overhead;
+
+
+ protected:
+
+ private:
+
+
+ void prepare_AR( int desired_blockSize, int n, int totalR, int columnsR);
+ void finalize_AR();
+
+
+ void prepare_Y(int desired_blockSize, int n, int totalY);
+ void finalize_Y();
+
+ void prepare_AL( int columns, int n);
+ void finalize_AL();
+
+ void prepare_B();
+ void finalize_B();
+
+
+ static void* async_io(void *ptr );
+
+ struct databel_fvi * load_databel_fvi( const char *path );
+ void free_databel_fvi( struct databel_fvi **fvi );
+ FILE * fgls_fopen( const char *path, const char *mode );
+
+ void * fgls_malloc_impl( const char* file, long line, size_t size );
+
+ type_fileh FHandler;
+ type_fileh* Fhandler;
+
+
+
+
+
+
+
+
+
+
+
+};
+
+#endif // AIOWRAPPER_H
Added: pkg/OmicABELnoMM/Algorithm.cpp
===================================================================
--- pkg/OmicABELnoMM/Algorithm.cpp (rev 0)
+++ pkg/OmicABELnoMM/Algorithm.cpp 2013-12-16 08:37:15 UTC (rev 1469)
@@ -0,0 +1,1714 @@
+#include "Algorithm.h"
+
+Algorithm::Algorithm()
+{
+ //ctor
+}
+
+Algorithm::~Algorithm()
+{
+ //dtor
+}
+
+
+void Algorithm::solve(struct Settings params, struct Outputs &out, int type)
+{
+ switch (type)
+ {
+ case FULL_NEQ:
+ fullNEQ(params,out);
+ break;
+ case P_NEQ:
+ partialNEQ(params,out);
+ break;
+ case P_NEQ_B_OPT:
+ partialNEQ_Blocked_STL(params,out);
+ break;
+ case FULL_QR:
+ fullQR(params,out);
+ break;
+ case P_QR:
+ partialQR(params,out);
+ break;
+ case P_QR_B_OPT:
+ partialQR_Blocked_Rtl(params,out);
+ break;
+ case P_NEQ_B_OPT_MD:
+ partialNEQ_Blocked_STL_MD(params,out);
+ break;
+
+ default:
+ break;
+ }
+
+}
+
+void Algorithm::extract_subMatrix(type_precision* source, type_precision* dest,int dim1_source, int dim2_source,
+ int dim1_ini,int dim1_end,int dim2_ini,int dim2_end)
+{
+
+ int i,j,idx=0;
+ int size, source_ini;
+ for(i = dim2_ini; i<dim2_end; i++)
+ {
+ j = dim1_ini;
+ source_ini = i*dim1_source+j;
+ size = dim1_end-dim1_ini;
+ memcpy( (type_precision*)&dest[idx], (type_precision*)&source[source_ini], size * sizeof(type_precision) );
+// for(j = dim1_ini; j<dim1_end; j++)
+// {
+// dest[idx] = source[i*dim1_source+j];
+// idx++;
+// }
+ idx += size;
+ }
+
+}
+
+void Algorithm::prepare_Bfinal(type_precision* bfinal, type_precision* top,type_precision* bot,int dim1_b, int dim2_b,int dim1_b_bot)
+{
+ //memcpy are faster version of the fors
+ int i,k,w,top_idx,bot_idx,max = dim1_b*dim2_b;
+ int size;
+ top_idx = 0;
+ bot_idx = 0;
+ for(k = 0; k < dim2_b; k++)
+ {
+ size = k*dim1_b+(dim1_b-dim1_b_bot)-(k*dim1_b);
+ memcpy( (type_precision*)&bfinal[k*dim1_b], (type_precision*)&top[top_idx], size * sizeof(type_precision) );
+// for(i = k*dim1_b; i < k*dim1_b+(dim1_b-dim1_b_bot); i++)
+// {
+// bfinal[i] = top[top_idx];
+// top_idx++;
+// }
+ top_idx += size;
+ i = k*dim1_b + size;
+ w=i;
+
+ size = w+dim1_b_bot - w;
+ memcpy( (type_precision*)&bfinal[w], (type_precision*)&bot[bot_idx], size * sizeof(type_precision) );
+// for(i = w; i < w+dim1_b_bot; i++)
+// {
+// bfinal[i] = bot[bot_idx];
+// bot_idx++;
+// }
+ bot_idx += size;
+ }
+
+}
+
+void Algorithm::prepare_QY(type_precision* qy, type_precision* top,type_precision* bot,int dim1_QY, int dim2_QY,int dim1_qy_bot,int bot_blocks )
+{
+
+ int i,k,w,top_idx,bot_idx,max = dim1_QY*dim2_QY;
+ top_idx = 0;
+ bot_idx = 0;
+ for(k = 0; k < dim2_QY; k++)
+ {
+ for(i = k*dim1_QY; i < (k+1)*dim1_QY-dim1_qy_bot; i++)
+ {
+ qy[i] = top[top_idx];
+ top_idx++;
+ }
+ w=i;
+
+ for(i = w; i < w+dim1_qy_bot; i++)
+ {
+ qy[i] = bot[bot_idx];
+ bot_idx++;
+ }
+ bot_idx+=(bot_blocks-1)*dim1_qy_bot;
+ }
+
+}
+
+type_precision* Algorithm::extract_R(type_precision* A,int dim1_A, int dim2_A)
+{
+ type_precision* R = (type_precision*)calloc(dim2_A*dim2_A,sizeof(type_precision));
+ int i,j;
+
+ int R_idx=0;
+
+ for(i = 0; i < dim2_A; i++)
+ {
+ for(j = 0; j <= i; j++)
+ {
+ R[R_idx] = A[j+i*dim1_A];
+ R_idx++;
+ }
+ R_idx = dim2_A*(i+1);
+ }
+ return R;
+}
+
+type_precision* Algorithm::prepare_R(type_precision* RL,int dim1_A, int dim2_AL,int dim2_AR)
+{
+ int R_dim = (dim2_AR+dim2_AL);
+ type_precision* R = new type_precision[R_dim*R_dim];
+
+ int i,j;
+
+ int RL_idx=0;
+ int R_idx=0;
+
+ for(i = 0; i < dim2_AL; i++)
+ {
+ for(j = 0; j <= i; j++)
+ {
+ RL_idx = i*dim1_A + j;
+ R_idx = i*R_dim + j;
+ R[R_idx] = RL[RL_idx];
+
+ }
+
+ }
+ return R;
+}
+
+void Algorithm::update_R(type_precision* R, type_precision* topRr, type_precision* botRr,int dim1, int dim2, int r)
+{
+ int i,j,w;
+ int max = dim1*dim2;
+ int rtr_idx=0;
+ int rbr_idx=0;
+ for( j = r; j > 0 ; j--)
+ {
+ for(i = max-dim1*j; i < max-dim1*j+dim2-r; i++)
+ {
+ R[i] = topRr[rtr_idx];
+ rtr_idx++;
+ }
+ w = i;
+ for(i = w; i < w+r; i++)
+ {
+ R[i] = botRr[rbr_idx];
+ rbr_idx++;
+ }
+
+ }
+
+}
+
+
+void Algorithm::build_S(type_precision* S,type_precision* Stl,type_precision* Str,type_precision* Sbr,int l,int r)
+{
+ int Sidx;
+ int p = l+r;
+ for(int i = 0; i < l; i++)
+ {
+ Sidx = i*p;
+ for(int j= 0; j <= i; j++)
+ {
+ S[Sidx] = Stl[j+i*l];
+ Sidx++;
+ }
+ }
+ for(int i = 0; i < r; i++)
+ {
+ Sidx = l*p+p*i;
+ for(int j= 0; j < l; j++)
+ {
+ S[Sidx] = Str[j+i*l];
+ Sidx++;
+ }
+ }
+
+ for(int i = 0; i < r; i++)
+ {
+ Sidx = l*p+l+p*i;
+ for(int j= 0; j <= i; j++)
+ {
+ S[Sidx] = Sbr[j+i*r];
+ Sidx++;
+ }
+ }
+
+
+
+}
+
+void Algorithm::check_result(type_precision* AL, type_precision* AR,int rowsA,int colsA, int rhs,int colsAR,
+ type_precision* y, type_precision* res)
+{
+ type_precision* A = (type_precision*)malloc(rowsA*colsA*sizeof(type_precision));
+
+ int i,ar_idx=0;
+ for(i = 0; i <rowsA*(colsA-colsAR) ; i++)
+ {
+ A[i] = AL[i];
+ }
+
+ for(i = rowsA*(colsA-colsAR); i <rowsA*colsA ; i++)
+ {
+ A[i] = AR[ar_idx];
+ ar_idx++;
+ }
+
+
+
+ type_precision* ynew = replicate_vec(y,rowsA*rhs);
+ type_precision* new_sol = (type_precision*)malloc(colsA*rhs*sizeof(type_precision));
+
+ lapack_int info = LAPACKE_sgels(STORAGE_TYPE,'N',rowsA,colsA,rhs,A,rowsA,ynew,rowsA);
+ assert(info == 0,"Error Check");
+
+
+
+
+
+ int index=0;
+ int index_new = 0;
+ for(i = 0; i < rhs; i++)
+ {
+ copy_vec( &ynew[index], &new_sol[index_new], colsA);
+ index += rowsA;
+ index_new += colsA;
+ }
+
+
+// if(PRINT)
+// printf("\n Btop=(Rtl\\(Ql'*Y))-(Rtl\\Rtr)*(Rbr\\(Qr'*Y)); \n [Btop ; Rbr\\Qr'*Y] - bcomputed \n");
+
+ cblas_saxpy(rhs*colsA, -1.0, res, 1,new_sol,1);
+ type_precision u_norm=cblas_snrm2(rhs*colsA,new_sol,1);
+ //
+ if(abs(u_norm) >= 0.0001 || isnan(u_norm))
+ {
+
+
+ fflush(stdout);
+ matlab_print_matrix("AL",rowsA,colsA-colsAR,AL);
+ matlab_print_matrix("AR",rowsA,colsAR,AR);
+ matlab_print_matrix("Y",rowsA,rhs,y);
+ printf("\nA = [AL AR]; [Q,R] = qr(A,0); rr=R\\(Q'*Y)\n");
+ matlab_print_matrix("bcomputed",colsA,rhs,res);
+ matlab_print_matrix("newsol",colsA,rhs,ynew);
+ printf("\n%%\tnrom: %0.2g", u_norm);
+ exit(1);
+ }
+ else
+ {
+ matlab_print_matrix("bcomputed",colsA,rhs,res);
+ matlab_print_matrix("newsol",colsA,rhs,ynew);
+ printf("\n%%\tnrom: %0.2g", u_norm);
+ }
+
+
+
+ //cout << "\t**************";
+
+
+ free(ynew);
+ free(new_sol);
+ free(A);
+
+
+
+
+
+
+}
+
+///////////////////////////////
+
+
+
+void Algorithm::partialNEQ_Blocked_STL_MD(struct Settings params, struct Outputs &out)
+{
+ int max_threads = params.threads;
+
+
+ srand (time(NULL));
+
+ blas_set_num_threads(max_threads);
+
+
+ type_precision *Ytemp;
+ lapack_int info,n,lda,ldy,l,r,k,p;
+
+ int i,j,w;
+ int m;
+
+ cputime_type start_tick,start_tick2, end_tick;
+
+ AIOfile.initialize(params);
+
+ n = params.n; l=params.l; r=params.r; p = l+r;
+
+ int y_amount = params.t;
+ int y_block_size = params.tb;//kk
+
+ int a_amount = params.m;
+ int a_block_size = params.mb;
+
+ int a_iters = (a_amount + a_block_size - 1) / a_block_size;
+
+ int y_iters = (y_amount + y_block_size - 1) / y_block_size;
+
+
+ lda = n; ldy = n;
+ k = p;
+
+ for(int j = 0; j < y_iters; j++)
+ {
+
+ if(y_iters >= 40 && (i%(y_iters/40))==0)
+ {
+ cout << "*" << flush;
+
+ }
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 1469
More information about the Genabel-commits
mailing list