67 if (
hit_Rank == 0 ) fprintf(stderr,
" \nUsage: lu <matrixSize> <bloqSize>\n\n");
76 int main(
int argc,
char** argv) {
87 blockSize = atoi(argv[2]);
96 HitTile_double matrix, vector;
101 int blockSizes[2] = { blockSize, blockSize };
102 int blockSizesB[2] = { 1, blockSize };
128 randomMat( matrixSelection, n, n, 0xCAFE, -50.0, 50.0 );
131 randomMat( vectorSelection, 1, n, 0xBECA, -50.0, 50.0 );
135 print( matrixSelection,
"A.dtxt" );
147 solveMat( matrixLayout, matrixSelection, vectorSelection );
151 print( matrixSelection,
"LU.dtxt" );
179 #define matrixBlockAt(m,i,j) hit_blockTileElemAtNS( m, double, 2, i, j )
180 #define matrixBufferAt(m,i,j) hit_blockTileBufferAtNS( m, double, 2, i, j )
181 #define matrixBlockAtArrayCoords(m,i,j) hit_blockTileElemAtArrayCoords( m, double, 2, i, j )
193 if (min>=max)
return;
199 int indBlock[2] = { 0, 0 };
200 int ind[2] = { 0, 0 };
202 for (ind[0] = 0; ind[0]<n1; ind[0]++ ) {
204 indBlock[0] = ind[0] / selection.childSize[0];
206 for (i=0; i<n2; i++)
drand48();
211 int numBlocks = (int)ceil( (
double)n2 / selection.childSize[1] );
212 for (indBlock[1] = 0; indBlock[1] < numBlocks; indBlock[1]++ ) {
215 if ( indBlock[1] != numBlocks-1 ) size = selection.childSize[1];
216 else size = ( n2 % selection.childSize[1] == 0 ) ? selection.childSize[1] : n2 % selection.childSize[1];
228 int localInd = ind[0] % selection.childSize[0];
229 for (ind[1] = 0; ind[1] <
size; ind[1]++)
263 void strsmuv(HitTile_double a, HitTile_double b) {
269 if (m==0 || n==0)
return;
271 for(k=m-1;k>=0;k--) {
273 if((temp=-hit_tileGetNoStride(a,2,k,i))!=0.0) {
280 temp=hit_tileGetNoStride(a,2,k,k);
283 hit_tileSetNoStride(b,2,hit_tileGetNoStride(b,2,j,k)/temp,j,k);
295 void strsmlv(HitTile_double a, HitTile_double b) {
302 if (m==0 || n==0)
return;
306 if((temp=-hit_tileGetNoStride(a,2,k,i))!=0.0) {
321 void strsml2( HitTile_double a, HitTile_double b ) {
327 if (m==0 || n==0)
return;
331 if ( hit_tileGetNoStride(b,2,j,k) != 0.0 ) {
333 temp=-hit_tileGetNoStride(b,2,j,k);
345 void dgemm ( HitTile_double a, HitTile_double b, HitTile_double c ) {
353 if (n==0 || m==0 || k==0)
return;
357 if((temp=-hit_tileGetNoStride(a,2,i,l))!=0.0) {
371 void dgemmv ( HitTile_double a, HitTile_double b, HitTile_double c ) {
379 if (n==0 || m==0 || k==0)
return;
383 if((temp=-hit_tileGetNoStride(a,2,i,l))!=0.0) {
402 int stage, i, ii, j, jj, k;
403 HitTile_double
block;
406 HitTile_double bufferBlock;
420 int processedBlocks[2] = { 0, 0 };
425 for( stage=0; stage<numStages; stage++ ) {
427 int stageRow = processedBlocks[0];
428 int stageColumn = processedBlocks[1];
433 HitCom combloqfil, combloqfilrcv;
435 combloqfilrcv =
hit_comBroadcastDimSelect( layout, 0,
hit_layDimOwner( layout, 0, stage ), &bufferStageRow, hit_shape2(hit_sigIndex(0), hit_sig( processedBlocks[1],
hit_tileDimCard(mat,1)-1,1)),
HIT_COM_TILECOORDS,
HIT_DOUBLE );
438 combloqfil =
hit_comBroadcastDimSelect( layout, 0,
HIT_COM_MYSELF, &mat, hit_shape2( hit_sigIndex(stageRow), hit_sig(processedBlocks[1],
hit_tileDimCard(mat,1)-1,1)),
HIT_COM_TILECOORDS,
HIT_DOUBLE );
441 HitCom combloqcol, combloqcolrcv;
443 combloqcolrcv =
hit_comBroadcastDimSelect( layout, 1,
hit_layDimOwner( layout, 1, stage ), &bufferStageCol, hit_shape2( hit_sig( processedBlocks[0],
hit_tileDimCard(mat,0)-1,1), hit_sigIndex(0)),
HIT_COM_TILECOORDS,
HIT_DOUBLE );
446 combloqcol =
hit_comBroadcastDimSelect( layout, 1,
HIT_COM_MYSELF, &mat, hit_shape2( hit_sig( processedBlocks[0],
hit_tileDimCard(mat,0)-1,1), hit_sigIndex(stageColumn)),
HIT_COM_TILECOORDS,
HIT_DOUBLE );
457 for( k=0; k < numBlockRows-1; k++ ) {
458 double diagonalElement = hit_tileGetNoStride( block, 2, k, k );
459 if ( diagonalElement == 0.0 ) fprintf(stderr,
"Warning: 0.0 in position %d,%d\n",
463 if ( k == numBlockRows-1 )
break;
466 for( ii=k+1; ii < numBlockRows; ii++ ) {
468 double elem = hit_tileGetNoStride(block,2,ii,k);
469 if ( elem != 0.0 ) hit_tileSetNoStride( block, 2, elem/diagonalElement, ii, k );
476 hit_shape2( hit_sig( k+1, numBlockRows-1, 1 ), hit_sigIndex(k) ),
482 for( ii = k+1; ii < numBlockRows; ii++ ) {
483 double elem = - hit_tileGetNoStride( block, 2, ii, k );
487 for( jj=k+1; jj < numBlockColumns; jj++ ) {
493 HitTile_double updateBlock =
matrixBlockAt( mat, stageRow, j );
564 for( k=0; k < numRowsBlock-1; k++ ){
567 HitCom comColBellowDiagonal =
hit_comBroadcastDimSelect( layout, 1, owner, &block,
hit_shape(2, hit_sig( k+1, numRowsBlock-1, 1), hit_sigIndex(0) ),
HIT_COM_TILECOORDS,
HIT_DOUBLE );
572 for( ii = k+1; ii < numRowsBlock; ii++ ) {
573 double elem = - hit_tileGetNoStride( block, 2, ii, 0 );
576 HitTile_double updateBlock =
matrixBlockAt( mat, stageRow, j );
649 HitTile_double bufferVectorBlock;
652 if (
hit_laySelfRanksDim( layout, 0 ) == 0 ) comReduceVectorBlock =
hit_comReduceDimSelect( layout, 1,
HIT_RANKS_NULL, &bufferVectorBlock,
HIT_SHAPE_WHOLE,
HIT_COM_TILECOORDS, &bufferVectorBlock,
HIT_SHAPE_WHOLE,
HIT_COM_TILECOORDS,
HIT_DOUBLE, suma );
662 HitTile_double blockA, blockB;
663 HitCom combloqfil,combloqfilrcv;
665 int processedBlocks[2] = { 0, 0 };
673 for( stage=0; stage < numStages; stage++ ) {
674 int stageRow = processedBlocks[0];
675 int stageColumn = processedBlocks[1];
681 combloqfilrcv =
hit_comRecvSelect( layout, hit_ranks2(
hit_layDimOwner( layout, 0, stage ),
hit_laySelfRanksDim( layout, 1 ) ), &bufferStageRow, hit_shape2( hit_sigIndex(0), hit_sig( 0, processedBlocks[1]-1, 1 ) ),
HIT_COM_TILECOORDS,
HIT_DOUBLE);
691 for( j = stageColumn-1; j >= 0; j-- ) {
701 for( j = stageColumn-1; j >= 0; j-- ) {
727 for( j=processedBlocks[1]-1; j >= 0; j-- ) {
734 if( processedBlocks[1] > 0 )
hit_comDo( &combloqfilrcv );
736 for( j=processedBlocks[1]-1; j >= 0; j-- )
749 if ( processedBlocks[1] > 0 ) {
750 combloqfil =
hit_comSendSelect( layout, hit_ranks2( 0,
hit_laySelfRanksDim( layout, 1 ) ), &mat, hit_shape2( hit_sigIndex( stageRow ), hit_sig( 0, processedBlocks[1]-1, 1) ),
HIT_COM_TILECOORDS,
HIT_DOUBLE);
766 for( stage = numStages-1; stage >= 0; stage-- ) {
767 int stageRow = processedBlocks[0];
768 int stageColumn = processedBlocks[1];
774 combloqfilrcv =
hit_comRecvSelect( layout, hit_ranks2(
hit_layDimOwner( layout, 0, stage),
hit_laySelfRanksDim( layout, 1 ) ), &bufferStageRow, hit_shape2( hit_sigIndex(0), hit_sig( processedBlocks[1]+1,
hit_tileDimCard(mat,1)-1, 1 ) ),
HIT_COM_TILECOORDS,
HIT_DOUBLE);
838 combloqfil =
hit_comSendSelect( layout, hit_ranks2( 0,
hit_laySelfRanksDim( layout, 1 ) ), &mat, hit_shape2( hit_sigIndex( stageRow ), hit_sig( processedBlocks[1]+1,
hit_tileDimCard(mat,1)-1, 1) ),
HIT_COM_TILECOORDS,
HIT_DOUBLE);
#define hit_shape(nd,...)
#define hit_tileHasArrayCoords(var, ndims,...)
#define hit_layShape(lay)
int hit_layDimOwner(HitLayout lay, int dim, int ind)
void strsmuv(HitTile_double a, HitTile_double b)
#define hit_tileNewType(baseType)
HitCom comReduceVectorBlock
void hit_comOpSumDouble(void *, void *, int *, HitType *)
void randomMat(HitBlockTile selection, int n1, int n2, int seed, double min, double max)
#define hit_comOpFree(operation)
#define hit_comSendSelect(lay, sendTo, tileP, selection, modeSelect, baseType)
#define hit_comUpdateDimSendTo(comm, dim, SendTo)
void hit_comFree(HitCom issue)
void hit_comDo(HitCom *issue)
#define hit_tileDimSig(var, dim)
#define hit_comUpdateSendTile(comm, sendTile)
#define hit_tileDomainAlloc(newVarP, baseType, numDims,...)
#define hit_comUpdateRecvTile(comm, recvTile)
#define hit_topology(name,...)
#define hit_sigIn(sig, ind)
#define hit_blockTileSelect(...)
void print(HitBlockTile selection, const char *fileName)
#define hit_laySelfRanksDim(lay, dim)
#define hit_clockPrintMax(c)
#define matrixBlockAtArrayCoords(m, i, j)
#define hit_tileDimCard(var, dim)
void factorizeMat(HitLayout, HitBlockTile)
#define hit_blockTileDimBlockSize(tile, dim)
HitCom hit_comBroadcastDimSelect(HitLayout lay, int dim, int root, const void *tile, HitShape selection, int modeSelect, HitType baseType)
void dgemm(HitTile_double a, HitTile_double b, HitTile_double c)
#define hit_tileTextFileWrite(var, file, coord, datatype, s1, s2)
void hit_topFree(HitTopology topo)
#define hit_comOp(function, operation)
HitCom comBroadcastDiagonalBlock
void hit_blockTileFree(HitBlockTile container)
void dgemmv(HitTile_double a, HitTile_double b, HitTile_double c)
#define hit_tileDimHasArrayCoord(var, dim, pos)
#define hit_tileTile2Array(var, dim, pos)
void hit_comInit(int *pargc, char **pargv[])
#define hit_layImActive(lay)
void solveMat(HitLayout, HitBlockTile, HitBlockTile)
#define hit_tileDomain(newVarP, baseType, numDims,...)
#define hit_tileElemAtNoStride(var, ndims,...)
void hit_layFree(HitLayout lay)
int main(int argc, char *argv[])
#define hit_blockTileAlloc(varP, baseType)
#define hit_comRecvSelect(lay, receiveFrom, tileP, selection, modeSelect, baseType)
#define matrixBufferAt(m, i, j)
#define matrixBlockAt(m, i, j)
#define hit_tileForDimDomain(tile, dim, index)
HitCom hit_comReduceDimSelect(HitLayout lay, int dim, HitRanks root, const void *tilePSend, HitShape selectionSend, int modeSelectSend, const void *tilePRecv, HitShape selectionRecv, int modeSelectRecv, HitType baseType, HitOp operation)
#define hit_comAllowDims(lay)
#define hit_layout(name, topo,...)
#define HIT_COM_TILECOORDS
void strsmlv(HitTile_double a, HitTile_double b)
#define hit_tileFree(var)
void hit_blockTileNew(HitBlockTile *container, void *originalVar, int blockSizes[HIT_MAXDIMS])
void strsml2(HitTile_double a, HitTile_double b)
#define hit_tileShape(var)
#define hit_clockReduce(lay, c)
#define hit_clockStart(c)