src/swsse2/swwozniak.c

   1 /******************************************************************\r
   2   Copyright 2006 by Michael Farrar.  All rights reserved.\r
   3   This program may not be sold or incorporated into a commercial product,\r
   4   in whole or in part, without written consent of Michael Farrar.  For \r
   5   further information regarding permission for use or reproduction, please \r
   6   contact: Michael Farrar at farrar.michael@gmail.com.\r
   7 *******************************************************************/\r
   8 \r
   9 /* Written by Michael Farrar, 2006.\r
  10    Please send bug reports and/or suggestions to farrar.michael@gmail.com.\r
  11 */\r
  12 \r
  13 /* Implementation of the Wozniak "vertical" vectorization\r
  14    strategy for Smith-Waterman comparison, Wozniak (1997) Comp.\r
  15    Appl. Biosci. 13:145-150\r
  16 */\r
  17 \r
  18 #include <stdio.h>\r
  19 #include <stdlib.h>\r
  20 #include <emmintrin.h>\r
  21 \r
  22 #include "swsse2.h"\r
  23 #include "swwozniak.h"\r
  24 \r
  25 #define MATRIX_ROW_SIZE 32\r
  26 #define MATRIX_SIZE     (MATRIX_ROW_SIZE * (ALPHA_SIZE + 1))\r
  27 \r
  28 typedef struct {\r
  29     char           *pbMatrix;\r
  30     short          *psMatrix;\r
  31     __m128i        *pvHStore;\r
  32     __m128i        *pvEStore;\r
  33     unsigned char  *pData;\r
  34     unsigned short  bias;\r
  35 } SwWozniakData;\r
  36 \r
  37 int\r
  38 swWozniakWord (unsigned char   *querySeq,\r
  39                int              queryLength,\r
  40                unsigned char   *dbSeq,\r
  41                int              dbLength,\r
  42                unsigned short   gapOpen,\r
  43                unsigned short   gapExtend,\r
  44                short           *pMatrix,\r
  45                __m128i         *pvHStore,\r
  46                __m128i         *pvEStore);\r
  47 \r
  48 int\r
  49 swWozniakByte (unsigned char   *querySeq,\r
  50                int              queryLength,\r
  51                unsigned char   *dbSeq,\r
  52                int              dbLength,\r
  53                unsigned short   gapOpen,\r
  54                unsigned short   gapExtend,\r
  55                char            *pMatrix,\r
  56                __m128i         *pvHStore,\r
  57                __m128i         *pvEStore,\r
  58                unsigned short   bias);\r
  59 \r
  60 void *\r
  61 swWozniakInit(unsigned char   *querySeq,\r
  62               int              queryLength,\r
  63               signed char     *matrix)\r
  64 {\r
  65     int i, j;\r
  66 \r
  67     int nCount;\r
  68 \r
  69     int lenQryByte;\r
  70     int lenQryShort;\r
  71 \r
  72     int bias;\r
  73 \r
  74     int weight;\r
  75 \r
  76     short *ps;\r
  77     char *pc;\r
  78 \r
  79     signed char *matrixRow;\r
  80 \r
  81     size_t aligned;\r
  82 \r
  83     SwWozniakData *pSwData;\r
  84  \r
  85     lenQryByte = (queryLength + 15) / 16 + 2;\r
  86     lenQryShort = (queryLength + 7) / 8 + 2;\r
  87 \r
  88     pSwData = (SwWozniakData *) malloc (sizeof (SwWozniakData));\r
  89     if (!pSwData) {\r
  90         fprintf (stderr, "Unable to allocate memory for SW data\n");\r
  91         exit (-1);\r
  92     }\r
  93 \r
  94     nCount = 64 +                             /* slack bytes */\r
  95              4 * (ALPHA_SIZE + 1) +           /* byte matrix */\r
  96              4 * (ALPHA_SIZE + 1) +           /* short matrix */\r
  97              ((queryLength + 16) * 2);        /* vH and vE */\r
  98 \r
  99 \r
 100     pSwData->pData = (unsigned char *) calloc (nCount, sizeof (__m128i));\r
 101     if (!pSwData->pData) {\r
 102         fprintf (stderr, "Unable to allocate memory for SW data buffers\n");\r
 103         exit (-1);\r
 104     }\r
 105 \r
 106     pSwData->pbMatrix = (char *) pSwData->pData;\r
 107     pSwData->psMatrix = (short *) (pSwData->pbMatrix + MATRIX_SIZE);\r
 108 \r
 109     /* since we might port this to another platform, lets align the data */\r
 110     /* to 16 byte boundries ourselves */\r
 111     aligned = (size_t) (pSwData->psMatrix + MATRIX_SIZE);\r
 112     aligned = (aligned + 15) & ~(0x0f);\r
 113 \r
 114     pSwData->pvHStore = (__m128i *) aligned;\r
 115     pSwData->pvEStore = pSwData->pvHStore + queryLength + 16;\r
 116 \r
 117     /* Find the bias to use in the substitution matrix */\r
 118     bias = 127;\r
 119     for (i = 0; i < ALPHA_SIZE * ALPHA_SIZE; i++) {\r
 120         if (matrix[i] < bias) {\r
 121             bias = matrix[i];\r
 122         }\r
 123     }\r
 124     if (bias > 0) {\r
 125         bias = 0;\r
 126     }\r
 127 \r
 128     pc = pSwData->pbMatrix;\r
 129     ps = pSwData->psMatrix;\r
 130 \r
 131     for (i = 0; i < ALPHA_SIZE; i++) {\r
 132         matrixRow = matrix + i * ALPHA_SIZE;\r
 133 \r
 134         for (j = 0; j < ALPHA_SIZE; j++) {\r
 135             weight = matrixRow[j];\r
 136             *pc++ = weight - bias;\r
 137             *ps++ = weight;\r
 138         }\r
 139 \r
 140         for ( ; j < MATRIX_ROW_SIZE; j++) {\r
 141             *pc++ = -bias;\r
 142             *ps++ = 0;\r
 143         }\r
 144     }\r
 145 \r
 146     /* add the weights for the NULL rows */\r
 147     for (j = 0; j < MATRIX_ROW_SIZE; j++) {\r
 148         *pc++ = -bias;\r
 149         *ps++ = 0;\r
 150     }\r
 151 \r
 152     pSwData->bias = (unsigned short) -bias;\r
 153 \r
 154     return pSwData;\r
 155 }\r
 156 \r
 157 void swWozniakScan (unsigned char   *querySeq,\r
 158                     int              queryLength,\r
 159                     FASTA_LIB       *dbLib,\r
 160                     void            *swData,\r
 161                     SEARCH_OPTIONS  *options,\r
 162                     SCORE_LIST      *scores)\r
 163 {\r
 164     int score;\r
 165 \r
 166     int threshold = options->threshold;\r
 167 \r
 168     unsigned char *dbSeq;\r
 169     int dbLen;\r
 170 \r
 171     int gapInit = -(options->gapInit + options->gapExt);\r
 172     int gapExt  = -options->gapExt;\r
 173 \r
 174     SwWozniakData *wozniakData = (SwWozniakData *) swData;\r
 175 \r
 176     dbSeq = nextSeq (dbLib, &dbLen);\r
 177     while (dbLen > 0) {\r
 178 \r
 179         score = swWozniakByte (querySeq, queryLength, \r
 180                                dbSeq, dbLen, \r
 181                                gapInit, gapExt, \r
 182                                wozniakData->pbMatrix,\r
 183                                wozniakData->pvHStore,\r
 184                                wozniakData->pvEStore,\r
 185                                wozniakData->bias);\r
 186 \r
 187         /* check if needs a run with higher precision */\r
 188         if (score >= 255) {\r
 189             score = swWozniakWord (querySeq, queryLength, \r
 190                                    dbSeq, dbLen, \r
 191                                    gapInit, gapExt, \r
 192                                    wozniakData->psMatrix,\r
 193                                    wozniakData->pvHStore,\r
 194                                    wozniakData->pvEStore);\r
 195         }\r
 196 \r
 197         if (score >= threshold) {\r
 198             int minScore = insertList (scores, score, seqName (dbLib));\r
 199             if (minScore >= threshold) {\r
 200                 threshold = minScore;\r
 201             }\r
 202         }\r
 203 \r
 204         dbSeq = nextSeq (dbLib, &dbLen);\r
 205     }\r
 206 }\r
 207 \r
 208 void\r
 209 swWozniakComplete(void *pSwData)\r
 210 {\r
 211     SwWozniakData *pWozniakData = (SwWozniakData *) pSwData;\r
 212 \r
 213     free (pWozniakData->pData);\r
 214     free (pWozniakData);\r
 215 }\r
 216 \r
 217 \r
 218 int\r
 219 swWozniakWord (unsigned char   *querySeq,\r
 220                int              queryLength,\r
 221                unsigned char   *dbSeq,\r
 222                int              dbLength,\r
 223                unsigned short   gapOpen,\r
 224                unsigned short   gapExtend,\r
 225                short           *pMatrix,\r
 226                __m128i         *pvHStore,\r
 227                __m128i         *pvEStore)\r
 228 {\r
 229     int      i, j, k, l, m;\r
 230     int      score;\r
 231 \r
 232     short   *pScore;\r
 233 \r
 234     __m128i  vE, vF, vH;\r
 235     __m128i  vEUp, vHUp1, vHUp2;\r
 236 \r
 237     __m128i  vMaxScore;\r
 238     __m128i  vGapOpen;\r
 239     __m128i  vGapExtend;\r
 240     __m128i  vScore;\r
 241 \r
 242     __m128i  vMin;\r
 243     __m128i  vMinimums;\r
 244     __m128i  vTemp;\r
 245 \r
 246     /* remove unreferenced warning */\r
 247     querySeq;\r
 248 \r
 249     /* Load gap opening penalty to all elements of a constant */\r
 250     vGapOpen = _mm_insert_epi16 (vGapOpen, gapOpen, 0);\r
 251     vGapOpen = _mm_shufflelo_epi16 (vGapOpen, 0);\r
 252     vGapOpen = _mm_shuffle_epi32 (vGapOpen, 0);\r
 253 \r
 254     /* Load gap extension penalty to all elements of a constant */\r
 255     vGapExtend = _mm_insert_epi16 (vGapExtend, gapExtend, 0);\r
 256     vGapExtend = _mm_shufflelo_epi16 (vGapExtend, 0);\r
 257     vGapExtend = _mm_shuffle_epi32 (vGapExtend, 0);\r
 258 \r
 259     /*  load vMaxScore with the zeros.  since we are using signed */\r
 260     /*  math, we will bias the maxscore to -32768 so we have the */\r
 261     /*  full range of the short. */\r
 262     vMaxScore = _mm_cmpeq_epi16 (vMaxScore, vMaxScore);\r
 263     vMaxScore = _mm_slli_epi16 (vMaxScore, 15);\r
 264 \r
 265     vMinimums = _mm_shuffle_epi32 (vMaxScore, 0);\r
 266 \r
 267     vMin = _mm_shuffle_epi32 (vMaxScore, 0);\r
 268     vMin = _mm_srli_si128 (vMin, 14);\r
 269 \r
 270     for (i = 0; i < queryLength + 8; i++)\r
 271     {\r
 272         _mm_store_si128 (pvEStore + i, vMaxScore);\r
 273         _mm_store_si128 (pvHStore + i, vMaxScore);\r
 274     }\r
 275 \r
 276     pScore = (short *) &vScore;\r
 277 \r
 278     for (i = 0; i < dbLength; i += 8, dbSeq += 8)\r
 279     {\r
 280         /* zero lots of stuff. */\r
 281         vE     = _mm_shuffle_epi32 (vMinimums, 0);\r
 282         vF     = _mm_shuffle_epi32 (vMinimums, 0);\r
 283         vH     = _mm_shuffle_epi32 (vMinimums, 0);\r
 284         vHUp2  = _mm_shuffle_epi32 (vMinimums, 0);\r
 285 \r
 286         vScore = _mm_xor_si128 (vScore, vScore);\r
 287 \r
 288         for (j = 0; j < 8; ++j)\r
 289         {\r
 290             for (k = 0; k <= j; ++k) {\r
 291                 int matrixOffset = *(dbSeq + k) * MATRIX_ROW_SIZE;\r
 292                 pScore[k] = *(pMatrix + matrixOffset + *(querySeq + j - k));\r
 293             }\r
 294             for ( ; k < 8; ++k) {\r
 295                 pScore[k] = 0;\r
 296             }\r
 297 \r
 298             /* load values of vE and vH from previous row (one unit up) */\r
 299             vEUp    = _mm_load_si128 (pvEStore + j);\r
 300             vHUp1   = _mm_load_si128 (pvHStore + j);\r
 301 \r
 302             /* shift into place so we have complete vE and vH vectors */\r
 303             /* that refer to the values one unit up from each cell */\r
 304             /* that we are currently working on. */\r
 305             vTemp   = _mm_slli_si128 (vE, 2);\r
 306             vEUp    = _mm_srli_si128 (vEUp, 14);\r
 307             vEUp    = _mm_or_si128 (vEUp, vTemp);\r
 308 \r
 309             vTemp   = _mm_slli_si128 (vH, 2);\r
 310             vHUp1   = _mm_srli_si128 (vHUp1, 14);\r
 311             vHUp1   = _mm_or_si128 (vHUp1, vTemp);\r
 312 \r
 313             /* do the dynamic programming */\r
 314 \r
 315             /* update vE value */\r
 316             vE    = _mm_subs_epi16 (vEUp, vGapExtend);\r
 317             vTemp = _mm_subs_epi16 (vHUp1, vGapOpen);\r
 318             vE    = _mm_max_epi16 (vE, vTemp);\r
 319 \r
 320             /* update vF value */\r
 321             vF    = _mm_subs_epi16 (vF, vGapExtend);\r
 322             vTemp = _mm_subs_epi16 (vH, vGapOpen);\r
 323             vF    = _mm_max_epi16  (vF, vTemp);\r
 324 \r
 325             /* add score to vH */\r
 326             vH = _mm_adds_epi16 (vHUp2, vScore);\r
 327 \r
 328             /* set vH to max of vH, vE, vF */\r
 329             vH = _mm_max_epi16 (vH, vE);\r
 330             vH = _mm_max_epi16 (vH, vF);\r
 331 \r
 332             /* Save value to use for next diagonal vH */\r
 333             vHUp2 = vHUp1;\r
 334 \r
 335             /* Update highest score encountered this far */\r
 336             vMaxScore = _mm_max_epi16 (vMaxScore, vH);\r
 337         }\r
 338 \r
 339         for (l = 0; j < queryLength; ++j, ++l)\r
 340         {\r
 341             for (k = 0; k < 8; ++k) {\r
 342                 int matrixOffset = *(dbSeq + k) * MATRIX_ROW_SIZE;\r
 343                 pScore[k] = *(pMatrix + matrixOffset + *(querySeq + j - k));\r
 344             }\r
 345 \r
 346             /* load values of vE and vH from previous row (one unit up) */\r
 347             vEUp    = _mm_load_si128 (pvEStore + j);\r
 348             vHUp1   = _mm_load_si128 (pvHStore + j);\r
 349 \r
 350             /* save old values of vE and vH to use on next row */\r
 351             _mm_store_si128 (pvEStore + l, vE);\r
 352             _mm_store_si128 (pvHStore + l, vH);\r
 353             \r
 354             /* shift into place so we have complete vE and vH vectors */\r
 355             /* that refer to the values one unit up from each cell */\r
 356             /* that we are currently working on. */\r
 357             vTemp = _mm_slli_si128 (vE, 2);\r
 358             vEUp  = _mm_srli_si128 (vEUp, 14);\r
 359             vEUp  = _mm_or_si128 (vEUp, vTemp);\r
 360 \r
 361             vTemp = _mm_slli_si128 (vH, 2);\r
 362             vHUp1 = _mm_srli_si128 (vHUp1, 14);\r
 363             vHUp1 = _mm_or_si128 (vHUp1, vTemp);\r
 364 \r
 365             /* do the dynamic programming */\r
 366 \r
 367             /* update vE value */\r
 368             vE    = _mm_subs_epi16 (vEUp, vGapExtend);\r
 369             vTemp = _mm_subs_epi16 (vHUp1, vGapOpen);\r
 370             vE    = _mm_max_epi16 (vE, vTemp);\r
 371 \r
 372             /* update vF value */\r
 373             vF    = _mm_subs_epi16 (vF, vGapExtend);\r
 374             vTemp = _mm_subs_epi16 (vH, vGapOpen);\r
 375             vF    = _mm_max_epi16 (vF, vTemp);\r
 376 \r
 377             /* add score to vH */\r
 378             vH = _mm_adds_epi16(vHUp2, vScore);\r
 379 \r
 380             /* set vH to max of vH, vE, vF */\r
 381             vH = _mm_max_epi16 (vH, vE);\r
 382             vH = _mm_max_epi16 (vH, vF);\r
 383 \r
 384             /* Save value to use for next diagonal vH */\r
 385             vHUp2 = vHUp1;\r
 386 \r
 387             /* Update highest score encountered this far */\r
 388             vMaxScore = _mm_max_epi16 (vMaxScore, vH);\r
 389         }\r
 390 \r
 391         for (m = 0 ; m < 7; ++j, ++l, ++m)\r
 392         {\r
 393             for (k = 0; k <= m; ++k) {\r
 394                 pScore[k] = 0;\r
 395             }\r
 396             for ( ; k < 8; ++k) {\r
 397                 int matrixOffset = *(dbSeq + k) * MATRIX_ROW_SIZE;\r
 398                 pScore[k] = *(pMatrix + matrixOffset + *(querySeq + j - k));\r
 399             }\r
 400 \r
 401             /* save old values of vE and vH to use on next row */\r
 402             _mm_store_si128 (pvEStore + l, vE);\r
 403             _mm_store_si128 (pvHStore + l, vH);\r
 404 \r
 405             /* v_score_load contains all zeros */\r
 406             vTemp = _mm_slli_si128 (vE, 2);\r
 407             vEUp  = _mm_or_si128 (vMin, vTemp);\r
 408             vTemp = _mm_slli_si128 (vH, 2);\r
 409             vHUp1 = _mm_or_si128 (vMin, vTemp);\r
 410 \r
 411             /* do the dynamic programming */\r
 412 \r
 413             /* update vE value */\r
 414             vE    = _mm_subs_epi16 (vEUp, vGapExtend);\r
 415             vTemp = _mm_subs_epi16 (vHUp1, vGapOpen);\r
 416             vE    = _mm_max_epi16 (vE, vTemp);\r
 417 \r
 418             /* update vF value */\r
 419             vF    = _mm_subs_epi16 (vF, vGapExtend);\r
 420             vTemp = _mm_subs_epi16 (vH, vGapOpen);\r
 421             vF    = _mm_max_epi16 (vF, vTemp);\r
 422 \r
 423             /* add score to vH */\r
 424             vH = _mm_adds_epi16 (vHUp2, vScore);\r
 425 \r
 426             /* set vH to max of vH, vE, vF */\r
 427             vH = _mm_max_epi16 (vH, vE);\r
 428             vH = _mm_max_epi16 (vH, vF);\r
 429 \r
 430             /* Save value to use for next diagonal vH */\r
 431             vHUp2 = vHUp1;\r
 432 \r
 433             /* Update highest score encountered this far */\r
 434             vMaxScore = _mm_max_epi16(vMaxScore,vH);\r
 435         }\r
 436 \r
 437         _mm_store_si128 (pvEStore + l, vE);\r
 438         _mm_store_si128 (pvHStore + l, vH);\r
 439     }\r
 440 \r
 441     /* find largest score in the vMaxScore vector */\r
 442     vTemp = _mm_srli_si128 (vMaxScore, 8);\r
 443     vMaxScore = _mm_max_epi16 (vMaxScore, vTemp);\r
 444     vTemp = _mm_srli_si128 (vMaxScore, 4);\r
 445     vMaxScore = _mm_max_epi16 (vMaxScore, vTemp);\r
 446     vTemp = _mm_srli_si128 (vMaxScore, 2);\r
 447     vMaxScore = _mm_max_epi16 (vMaxScore, vTemp);\r
 448 \r
 449     /* store in temporary variable */\r
 450     score = (short) _mm_extract_epi16 (vMaxScore, 0);\r
 451 \r
 452     /* return largest score */\r
 453     return score + SHORT_BIAS;\r
 454 }\r
 455 \r
 456 \r
 457 \r
 458 \r
 459 int\r
 460 swWozniakByte (unsigned char   *querySeq,\r
 461                int              queryLength,\r
 462                unsigned char   *dbSeq,\r
 463                int              dbLength,\r
 464                unsigned short   gapOpen,\r
 465                unsigned short   gapExtend,\r
 466                char            *pMatrix,\r
 467                __m128i         *pvHStore,\r
 468                __m128i         *pvEStore,\r
 469                unsigned short   bias)\r
 470 {\r
 471     int      i, j, k, l, m;\r
 472     int      score;\r
 473 \r
 474     int      dup;\r
 475 \r
 476     char    *pScore;\r
 477 \r
 478     __m128i  vE, vF, vH;\r
 479     __m128i  vEUp, vHUp1, vHUp2;\r
 480 \r
 481     __m128i  vMaxScore;\r
 482     __m128i  vBias;\r
 483     __m128i  vGapOpen;\r
 484     __m128i  vGapExtend;\r
 485     __m128i  vScore;\r
 486 \r
 487     __m128i  vTemp;\r
 488 \r
 489     /* remove unreferenced warning */\r
 490     querySeq;\r
 491 \r
 492     /* Load the bias to all elements of a constant */\r
 493     dup   = (bias << 8) | (bias & 0x00ff);\r
 494     vBias = _mm_insert_epi16 (vBias, dup, 0);\r
 495     vBias = _mm_shufflelo_epi16 (vBias, 0);\r
 496     vBias = _mm_shuffle_epi32 (vBias, 0);\r
 497 \r
 498     /* Load gap opening penalty to all elements of a constant */\r
 499     dup      = (gapOpen << 8) | (gapOpen & 0x00ff);\r
 500     vGapOpen = _mm_insert_epi16 (vGapOpen, dup, 0);\r
 501     vGapOpen = _mm_shufflelo_epi16 (vGapOpen, 0);\r
 502     vGapOpen = _mm_shuffle_epi32 (vGapOpen, 0);\r
 503 \r
 504     /* Load gap extension penalty to all elements of a constant */\r
 505     dup        = (gapExtend << 8) | (gapExtend & 0x00ff);\r
 506     vGapExtend = _mm_insert_epi16 (vGapExtend, dup, 0);\r
 507     vGapExtend = _mm_shufflelo_epi16 (vGapExtend, 0);\r
 508     vGapExtend = _mm_shuffle_epi32 (vGapExtend, 0);\r
 509 \r
 510     vScore = _mm_xor_si128 (vScore, vScore);\r
 511     vMaxScore = _mm_xor_si128 (vMaxScore, vMaxScore);\r
 512 \r
 513     for (i = 0; i < queryLength + 16; i++)\r
 514     {\r
 515         _mm_store_si128 (pvEStore + i, vMaxScore);\r
 516         _mm_store_si128 (pvHStore + i, vMaxScore);\r
 517     }\r
 518 \r
 519     pScore = (char *) &vScore;\r
 520 \r
 521     for (i = 0; i < dbLength; i += 16, dbSeq += 16)\r
 522     {\r
 523         // zero lots of stuff.\r
 524         vE     = _mm_xor_si128 (vE, vE);\r
 525         vF     = _mm_xor_si128 (vF, vF);\r
 526         vH     = _mm_xor_si128 (vH, vH);\r
 527         vHUp2  = _mm_xor_si128 (vHUp2, vHUp2);\r
 528 \r
 529         vScore = _mm_xor_si128 (vScore, vScore);\r
 530 \r
 531         for (j = 0; j < 16; ++j)\r
 532         {\r
 533             for (k = 0; k <= j; ++k) {\r
 534                 int matrixOffset = *(dbSeq + k) * MATRIX_ROW_SIZE;\r
 535                 pScore[k] = *(pMatrix + matrixOffset + *(querySeq + j - k));\r
 536             }\r
 537             for ( ; k < 16; ++k) {\r
 538                 pScore[k] = (char) bias;\r
 539             }\r
 540 \r
 541             // load values of vE and vH from previous row (one unit up)\r
 542             vEUp    = _mm_load_si128 (pvEStore + j);\r
 543             vHUp1   = _mm_load_si128 (pvHStore + j);\r
 544 \r
 545             // shift into place so we have complete vE and vH vectors\r
 546             // that refer to the values one unit up from each cell\r
 547             // that we are currently working on.\r
 548             vTemp   = _mm_slli_si128 (vE, 1);\r
 549             vEUp    = _mm_srli_si128 (vEUp, 15);\r
 550             vEUp    = _mm_or_si128 (vEUp, vTemp);\r
 551 \r
 552             vTemp   = _mm_slli_si128 (vH, 1);\r
 553             vHUp1   = _mm_srli_si128 (vHUp1, 15);\r
 554             vHUp1   = _mm_or_si128 (vHUp1, vTemp);\r
 555 \r
 556             // do the dynamic programming\r
 557 \r
 558             // update vE value\r
 559             vE    = _mm_subs_epu8 (vEUp, vGapExtend);\r
 560             vTemp = _mm_subs_epu8 (vHUp1, vGapOpen);\r
 561             vE    = _mm_max_epu8 (vE, vTemp);\r
 562 \r
 563             // update vF value\r
 564             vF    = _mm_subs_epu8 (vF, vGapExtend);\r
 565             vTemp = _mm_subs_epu8 (vH, vGapOpen);\r
 566             vF    = _mm_max_epu8  (vF, vTemp);\r
 567 \r
 568             // add score to vH\r
 569             vH = _mm_adds_epu8 (vHUp2, *((__m128i *) pScore));\r
 570             vH = _mm_subs_epu8 (vH, vBias);\r
 571 \r
 572             // set vH to max of vH, vE, vF\r
 573             vH = _mm_max_epu8 (vH, vE);\r
 574             vH = _mm_max_epu8 (vH, vF);\r
 575 \r
 576             // Save value to use for next diagonal vH\r
 577             vHUp2 = vHUp1;\r
 578 \r
 579             // Update highest score encountered this far\r
 580             vMaxScore = _mm_max_epu8 (vMaxScore, vH);\r
 581         }\r
 582 \r
 583         for (l = 0; j < queryLength; ++j, ++l)\r
 584         {\r
 585             for (k = 0; k < 16; ++k) {\r
 586                 int matrixOffset = *(dbSeq + k) * MATRIX_ROW_SIZE;\r
 587                 pScore[k] = *(pMatrix + matrixOffset + *(querySeq + j - k));\r
 588             }\r
 589 \r
 590             // load values of vE and vH from previous row (one unit up)\r
 591             vEUp    = _mm_load_si128 (pvEStore + j);\r
 592             vHUp1   = _mm_load_si128 (pvHStore + j);\r
 593 \r
 594             // save old values of vE and vH to use on next row\r
 595             _mm_store_si128 (pvEStore + l, vE);\r
 596             _mm_store_si128 (pvHStore + l, vH);\r
 597             \r
 598             // shift into place so we have complete vE and vH vectors\r
 599             // that refer to the values one unit up from each cell\r
 600             // that we are currently working on.\r
 601             vTemp = _mm_slli_si128 (vE, 1);\r
 602             vEUp  = _mm_srli_si128 (vEUp, 15);\r
 603             vEUp  = _mm_or_si128 (vEUp, vTemp);\r
 604 \r
 605             vTemp = _mm_slli_si128 (vH, 1);\r
 606             vHUp1 = _mm_srli_si128 (vHUp1, 15);\r
 607             vHUp1 = _mm_or_si128 (vHUp1, vTemp);\r
 608 \r
 609             // do the dynamic programming\r
 610 \r
 611             // update vE value\r
 612             vE    = _mm_subs_epu8 (vEUp, vGapExtend);\r
 613             vTemp = _mm_subs_epu8 (vHUp1, vGapOpen);\r
 614             vE    = _mm_max_epu8 (vE, vTemp);\r
 615 \r
 616             // update vF value\r
 617             vF    = _mm_subs_epu8 (vF, vGapExtend);\r
 618             vTemp = _mm_subs_epu8 (vH, vGapOpen);\r
 619             vF    = _mm_max_epu8 (vF, vTemp);\r
 620 \r
 621             // add score to vH\r
 622             vH = _mm_adds_epu8(vHUp2, vScore);\r
 623             vH = _mm_subs_epu8 (vH, vBias);\r
 624 \r
 625             // set vH to max of vH, vE, vF\r
 626             vH = _mm_max_epu8 (vH, vE);\r
 627             vH = _mm_max_epu8 (vH, vF);\r
 628 \r
 629             // Save value to use for next diagonal vH\r
 630             vHUp2 = vHUp1;\r
 631 \r
 632             // Update highest score encountered this far\r
 633             vMaxScore = _mm_max_epu8 (vMaxScore, vH);\r
 634         }\r
 635 \r
 636         for (m = 0 ; m < 15; ++j, ++l, ++m)\r
 637         {\r
 638             for (k = 0; k <= m; ++k) {\r
 639                 pScore[k] = (char) bias;\r
 640             }\r
 641             for ( ; k < 16; ++k) {\r
 642                 int matrixOffset = *(dbSeq + k) * MATRIX_ROW_SIZE;\r
 643                 pScore[k] = *(pMatrix + matrixOffset + *(querySeq + j - k));\r
 644             }\r
 645 \r
 646             // save old values of vE and vH to use on next row\r
 647             _mm_store_si128 (pvEStore + l, vE);\r
 648             _mm_store_si128 (pvHStore + l, vH);\r
 649 \r
 650             // v_score_load contains all zeros\r
 651             vEUp  = _mm_slli_si128 (vE, 1);\r
 652             vHUp1 = _mm_slli_si128 (vH, 1);\r
 653 \r
 654             // do the dynamic programming\r
 655 \r
 656             // update vE value\r
 657             vE    = _mm_subs_epu8 (vEUp, vGapExtend);\r
 658             vTemp = _mm_subs_epu8 (vHUp1, vGapOpen);\r
 659             vE    = _mm_max_epu8 (vE, vTemp);\r
 660 \r
 661             // update vF value\r
 662             vF    = _mm_subs_epu8 (vF, vGapExtend);\r
 663             vTemp = _mm_subs_epu8 (vH, vGapOpen);\r
 664             vF    = _mm_max_epu8 (vF, vTemp);\r
 665 \r
 666             // add score to vH\r
 667             vH = _mm_adds_epu8 (vHUp2, vScore);\r
 668             vH = _mm_subs_epu8 (vH, vBias);\r
 669 \r
 670             // set vH to max of vH, vE, vF\r
 671             vH = _mm_max_epu8 (vH, vE);\r
 672             vH = _mm_max_epu8 (vH, vF);\r
 673 \r
 674             // Save value to use for next diagonal vH\r
 675             vHUp2 = vHUp1;\r
 676 \r
 677             // Update highest score encountered this far\r
 678             vMaxScore = _mm_max_epu8(vMaxScore,vH);\r
 679         }\r
 680 \r
 681         _mm_store_si128 (pvEStore + l, vE);\r
 682         _mm_store_si128 (pvHStore + l, vH);\r
 683     }\r
 684 \r
 685     // find largest score in the vMaxScore vector\r
 686     vTemp = _mm_srli_si128 (vMaxScore, 8);\r
 687     vMaxScore = _mm_max_epu8 (vMaxScore, vTemp);\r
 688     vTemp = _mm_srli_si128 (vMaxScore, 4);\r
 689     vMaxScore = _mm_max_epu8 (vMaxScore, vTemp);\r
 690     vTemp = _mm_srli_si128 (vMaxScore, 2);\r
 691     vMaxScore = _mm_max_epu8 (vMaxScore, vTemp);\r
 692     vTemp = _mm_srli_si128 (vMaxScore, 1);\r
 693     vMaxScore = _mm_max_epu8 (vMaxScore, vTemp);\r
 694 \r
 695     // store in temporary variable\r
 696     score = (short) _mm_extract_epi16 (vMaxScore, 0);\r
 697     score = score & 0x00ff;\r
 698 \r
 699     //  check if we might have overflowed\r
 700     if (score + bias >= 255)\r
 701     {\r
 702         score = 255;\r
 703     }\r
 704 \r
 705 \r
 706     // return largest score\r
 707     return score;\r
 708 }\r