00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include "swps3.h"
00031 #include "matrix.h"
00032 #include "fasta.h"
00033 #include "DynProgr_scalar.h"
00034 #if defined(__SSE2__)
00035 #include "DynProgr_sse_byte.h"
00036 #include "DynProgr_sse_short.h"
00037 #endif
00038 #if defined(__ALTIVEC__)
00039 #include "DynProgr_altivec.h"
00040 #endif
00041 #if defined(__PS3)
00042 #include "DynProgr_PPU.h"
00043 #endif
00044 #include <stdio.h>
00045 #include <stdlib.h>
00046 #include <sys/types.h>
00047 #include <sys/wait.h>
00048 #include <sys/select.h>
00049 #include <unistd.h>
00050 #include <string.h>
00051 #include <assert.h>
00052 #include <float.h>
00053
00054 static ssize_t write_all(int fd, const void *data, size_t len) {
00055 size_t sent = 0;
00056 while(sent < len) {
00057 ssize_t res;
00058 res = write(fd,(const int8_t*)data+sent,len-sent);
00059 switch(res) {
00060 case 0:
00061 return sent;
00062 case -1:
00063 return -1;
00064 default:
00065 sent += res;
00066 }
00067 }
00068 return sent;
00069 }
00070
00071 static ssize_t read_all(int fd, void *buf, size_t len) {
00072 size_t recv = 0;
00073 while(recv < len) {
00074 ssize_t res;
00075 res = read(fd,(int8_t*)buf+recv,len-recv);
00076 switch(res) {
00077 case 0:
00078 return recv;
00079 case -1:
00080 return -1;
00081 default:
00082 recv += res;
00083 }
00084 }
00085 return recv;
00086 }
00087
00088 int main( int argc, char * argv[] ){
00089 char * matrixFile = NULL, * queryFile = NULL, * dbFile = NULL;
00090 int i, queryLen;
00091 char * query;
00092 #if defined(__SSE2__)
00093 SWType type = SSE2;
00094 #elif defined(__PS3)
00095 SWType type = PS3;
00096 #elif defined(__ALTIVEC__)
00097 SWType type = ALTIVEC;
00098 #else
00099 SWType type = SCALAR;
00100 #endif
00101 Options options = {-12,-2,DBL_MAX};
00102 int qCount=0, dCount=0, qResidues=0, dResidues=0;
00103 #ifdef HAVE_SYSCONF_NPROCESSORS
00104 int threads = sysconf(_SC_NPROCESSORS_ONLN);
00105 #else
00106 #if defined(__PS3)
00107 int threads = 6;
00108 #else
00109 int threads = 1;
00110 #endif
00111 #endif
00112 SBMatrix matrix;
00113 FastaLib * queryLib;
00114 for ( i=1; i<argc; i++ ){
00115 if (argv[i][0]=='-'){
00116 switch( argv[i][1] ){
00117 case 'h':
00118 matrixFile = NULL;
00119 i = argc; break;
00120 case 's':
00121 type = SCALAR;
00122 break;
00123 case 't':
00124 options.threshold = atoi( argv[++i] );
00125 break;
00126 case 'i':
00127 options.gapOpen = atoi( argv[++i] );
00128 break;
00129 case 'e':
00130 options.gapExt = atoi( argv[++i] );
00131 break;
00132 case 'j':
00133 threads = atoi( argv[++i] );
00134 break;
00135 default:
00136 matrixFile = NULL;
00137 i = argc; break;
00138 }
00139 }else{
00140 if (matrixFile == NULL)
00141 matrixFile = argv[i];
00142 else if (queryFile == NULL)
00143 queryFile = argv[i];
00144 else if (dbFile == NULL)
00145 dbFile = argv[i];
00146 else{
00147 matrixFile = NULL;
00148 i = argc; break;
00149 }
00150 }
00151 }
00152 if ( matrixFile == NULL || queryFile == NULL || dbFile == NULL ){
00153 printf( "Usage: %s [-h] [-s] [-j num] [-i num] [-e num] [-t num] matrix query db\n", argv[0] );
00154 return 0;
00155 }
00156 matrix = swps3_readSBMatrix( matrixFile );
00157 queryLib = swps3_openLib( queryFile );
00158 while ( (query=swps3_readNextSequence( queryLib, &queryLen )) ){
00159 double score = 0;
00160 int *pipes_read;
00161 int *pipes_write;
00162 char **seq_names;
00163 pid_t *children;
00164 int child_id = -1;
00165 int childpipe_read = -1;
00166 int childpipe_write = -1;
00167 FastaLib * dbLib = swps3_openLib( dbFile );
00168 #if defined(__SSE2__)
00169 ProfileByte * profileByte = swps3_createProfileByteSSE( query, queryLen, matrix );
00170 ProfileShort * profileShort = swps3_createProfileShortSSE( query, queryLen, matrix );
00171 #endif
00172 #if defined(__PS3)
00173 SPEProfile * profileByte = swps3_createProfileBytePPU(query, queryLen, matrix, MAX_SEQ_LENGTH);
00174 SPEProfile * profileShort = swps3_createProfileShortPPU(query, queryLen, matrix, MAX_SEQ_LENGTH);
00175 SPEProfile * profileFloat = swps3_createProfileFloatPPU(query, queryLen, matrix, MAX_SEQ_LENGTH);
00176
00177 int current_profile_is_byte = 0;
00178
00179 #endif
00180 #if defined(__ALTIVEC__)
00181 void *profileByteAltivec = swps3_createProfileByteAltivec(query, queryLen, matrix);
00182 void *profileShortAltivec = swps3_createProfileShortAltivec(query, queryLen, matrix);
00183 void *profileFloatAltivec = swps3_createProfileFloatAltivec(query, queryLen, matrix);
00184 #endif
00185 pipes_read = malloc(threads*sizeof(*pipes_read));
00186 pipes_write = malloc(threads*sizeof(*pipes_write));
00187 children = malloc(threads*sizeof(*children));
00188 seq_names = malloc(threads*sizeof(*seq_names));
00189 for(i=0;i<threads;++i) {
00190 pipes_read[i]=-1;
00191 pipes_write[i]=-1;
00192 children[i]=-1;
00193 seq_names[i]=malloc((MAX_SEQ_NAME_LENGTH+1)*sizeof(char));
00194 seq_names[i][MAX_SEQ_NAME_LENGTH]='\0';
00195 }
00196
00197 qCount++; qResidues+=queryLen;
00198 dCount=dResidues=0;
00199
00200 if(threads>1) {
00201 for(i=0; i<threads; ++i) {
00202 int fds[2];
00203 int res;
00204 char *db;
00205 int dbLen;
00206
00207 db=swps3_readNextSequence( dbLib, &dbLen);
00208 strncpy(seq_names[i],swps3_getSequenceName(dbLib),MAX_SEQ_NAME_LENGTH);
00209
00210 if(db == NULL) break;
00211
00212 dResidues+=dbLen;
00213
00214 res = pipe(fds);
00215 if(res < 0) {
00216 perror("error during pipe()");
00217 exit(1);
00218 }
00219 pipes_read[i] = fds[0];
00220 childpipe_write = fds[1];
00221
00222 res = pipe(fds);
00223 if(res < 0) {
00224 perror("error during pipe()");
00225 exit(1);
00226 }
00227 pipes_write[i] = fds[1];
00228 childpipe_read = fds[0];
00229
00230 children[i] = fork();
00231 if(children[i] < 0) {
00232 perror("error during fork()");
00233 exit(1);
00234 } else if(children[i] == 0) {
00235 int j;
00236 for(j=0; j<=i; ++j) {
00237 close(pipes_read[j]);
00238 close(pipes_write[j]);
00239 pipes_read[j] = -1;
00240 pipes_write[j] = -1;
00241 }
00242 child_id = i;
00243
00244 break;
00245 } else {
00246 ssize_t sres;
00247
00248 close(childpipe_read);
00249 close(childpipe_write);
00250 childpipe_read = -1;
00251 childpipe_write = -1;
00252
00253 sres = write_all(pipes_write[i],&dbLen,sizeof(int));
00254 if(sres != sizeof(int)) {
00255 perror("error during write()");
00256 exit(1);
00257 }
00258 sres = write_all(pipes_write[i],db,dbLen);
00259 if(sres != dbLen) {
00260 perror("error during write()");
00261 exit(1);
00262 }
00263 }
00264 }
00265 } else {
00266 childpipe_read = -1;
00267 childpipe_write = -1;
00268 child_id = 0;
00269 }
00270
00271 do {
00272 int dbLen;
00273 char * db;
00274 char * dbName;
00275
00276 if(childpipe_read <= 0) {
00277 db = swps3_readNextSequence( dbLib, &dbLen);
00278 dbName = swps3_getSequenceName(dbLib);
00279 } else {
00280 static char buffer[MAX_SEQ_LENGTH] __ALIGNED__;
00281 ssize_t res;
00282
00283 db = buffer;
00284
00285 res = read_all(childpipe_read,&dbLen,sizeof(int));
00286 if(res != sizeof(int)) {
00287 if(res == 0) exit(0);
00288 perror("error during read()");
00289 exit(1);
00290 }
00291 res = read_all(childpipe_read,db,dbLen);
00292 if(res != dbLen) {
00293 perror("error during read()");
00294 exit(1);
00295 }
00296 dbName = seq_names[child_id];
00297 }
00298
00299 if(child_id == -1) {
00300 fd_set readfds;
00301 int max_fd = -1;
00302 int i;
00303 int res;
00304
00305 FD_ZERO(&readfds);
00306
00307 for(i=0; i<threads; ++i) {
00308 if(pipes_read[i] > 0) FD_SET(pipes_read[i],&readfds);
00309 if(pipes_read[i]>max_fd) max_fd = pipes_read[i];
00310 }
00311
00312
00313 if(max_fd == -1) break;
00314
00315 res = select(max_fd+1, &readfds, NULL, NULL, NULL);
00316 if(res <= 0) {
00317 perror("error during select()");
00318 exit(1);
00319 }
00320 for(i=0; i<threads; ++i) {
00321 if(pipes_read[i] > 0 && FD_ISSET(pipes_read[i],&readfds)) {
00322 static char tmpbuff[MAX_SEQ_NAME_LENGTH+1];
00323 char *newName;
00324 res = read_all(pipes_read[i],&score,sizeof(score));
00325 if(res != sizeof(score)) {
00326 perror("error during read()");
00327 exit(1);
00328 }
00329
00330 newName = dbName;
00331 strncpy(tmpbuff,seq_names[i],MAX_SEQ_NAME_LENGTH);
00332 dbName = tmpbuff;
00333
00334 if(db) {
00335 strncpy(seq_names[i],newName,MAX_SEQ_NAME_LENGTH);
00336
00337 res = write_all(pipes_write[i],&dbLen,sizeof(int));
00338 if(res != sizeof(int)) {
00339 perror("error during write()");
00340 exit(1);
00341 }
00342 res = write_all(pipes_write[i],db,dbLen);
00343 if(res != dbLen) {
00344 perror("error during write()");
00345 exit(1);
00346 }
00347 } else {
00348 close(pipes_write[i]);
00349 close(pipes_read[i]);
00350 pipes_write[i] = -1;
00351 pipes_read[i] = -1;
00352 free(seq_names[i]);
00353 seq_names[i] = NULL;
00354 waitpid(children[i],NULL,0);
00355 children[i] = -1;
00356 }
00357 break;
00358 }
00359 }
00360 } else {
00361 if(db == NULL) break;
00362
00363 #ifdef DEBUG
00364 for(i=0; i<queryLen; ++i) printf("\t%c",query[i]);
00365 printf("\n");
00366 #endif
00367
00368 #if defined(__SSE2__)
00369 if(type == SSE2) {
00370 if( (score = swps3_alignmentByteSSE( profileByte, db, dbLen, &options )) >= DBL_MAX ) {
00371 score = swps3_alignmentShortSSE( profileShort, db, dbLen, &options );
00372 assert(score >= 250 && "score too low");
00373 }
00374 }
00375 #endif
00376 #if defined(__PS3)
00377 if(type == PS3) {
00378 #if defined(__ALTIVEC__)
00379 if(child_id == 6) {
00380 #if 0
00381 score = swps3_dynProgrFloatAltivec(db, dbLen, profileFloatAltivec, &options);
00382 #else
00383 score = swps3_dynProgrByteAltivec(db, dbLen, profileByteAltivec, &options);
00384 if(score >= DBL_MAX)
00385 score = swps3_dynProgrShortAltivec(db, dbLen, profileShortAltivec, &options);
00386 #endif
00387 } else {
00388 #endif
00389 #if 0
00390 loadProfileFloat(profileFloat, MAX_SEQ_LENGTH, &options);
00391 score = alignmentProfileSPE( db, dbLen );
00392 #elif 1
00393 if(!current_profile_is_byte) swps3_loadProfileByte(profileByte, MAX_SEQ_LENGTH, &options);
00394 current_profile_is_byte = 1;
00395 score = swps3_alignmentProfileSPE( db, dbLen );
00396 if( score >= DBL_MAX ) {
00397 if(current_profile_is_byte) swps3_loadProfileShort(profileShort, MAX_SEQ_LENGTH, &options);
00398 current_profile_is_byte = 0;
00399 score = swps3_alignmentProfileSPE( db, dbLen );
00400 }
00401 #elif 1
00402 if(!current_profile_is_byte) swps3_loadProfileShort(profileShort, MAX_SEQ_LENGTH, &options);
00403 current_profile_is_byte = 1;
00404 score = swps3_alignmentProfileSPE( db, dbLen );
00405 #else
00406 score = swps3_alignmentByteSPE( matrix, query, queryLen, db, dbLen, &options );
00407 if( score >= DBL_MAX ) {
00408 score = swps3_alignmentShortSPE( matrix, query, queryLen, db, dbLen, &options );
00409 }
00410 #endif
00411 #if defined(__ALTIVEC__)
00412 }
00413 #endif
00414 }
00415 #endif
00416 #if defined(__ALTIVEC__)
00417 if(type == ALTIVEC) {
00418 #if 0
00419 score = swps3_dynProgrFloatAltivec(db, dbLen, profileFloatAltivec, &options);
00420 #else
00421 score = swps3_dynProgrByteAltivec(db, dbLen, profileByteAltivec, &options);
00422 if(score >= DBL_MAX)
00423 score = swps3_dynProgrShortAltivec(db, dbLen, profileShortAltivec, &options);
00424 }
00425 #endif
00426 #endif
00427 if(type == SCALAR) {
00428 double dmatrix[MATRIX_DIM*MATRIX_DIM];
00429 for(i=0;i<MATRIX_DIM*MATRIX_DIM;++i) dmatrix[i]=matrix[i];
00430 score = swps3_alignScalar( dmatrix, query, queryLen, db, dbLen, &options);
00431 }
00432 }
00433
00434 if(childpipe_write > 0) {
00435 ssize_t res;
00436
00437 res = write_all(childpipe_write,&score,sizeof(score));
00438 if(res != sizeof(score)) {
00439 perror("error during write()");
00440 exit(1);
00441 }
00442 } else {
00443 if(score >= options.threshold) {
00444 printf(">threshold\t%s\n",dbName);
00445 } else {
00446 printf("%g\t%s\n",score,dbName);
00447 }
00448 }
00449
00450 dCount++; dResidues+=dbLen;
00451 } while(1);
00452
00453 free(pipes_read);
00454 free(pipes_write);
00455 free(children);
00456 free(seq_names);
00457
00458 #if defined(__SSE2__)
00459 swps3_freeProfileByteSSE( profileByte );
00460 swps3_freeProfileShortSSE( profileShort );
00461 #endif
00462 #if defined(__PS3)
00463 swps3_freeProfilePPU( profileByte );
00464 swps3_freeProfilePPU( profileShort );
00465 swps3_freeProfilePPU( profileFloat );
00466 #endif
00467 swps3_closeLib( dbLib );
00468 }
00469 fprintf(stderr,"%d[%d] x %d[%d]\n", qCount, qResidues, dCount, dResidues );
00470
00471 swps3_closeLib( queryLib );
00472 return 0;
00473 }
