+++ /dev/null
-/******************************************************************
- Copyright 2006 by Michael Farrar. All rights reserved.
- This program may not be sold or incorporated into a commercial product,
- in whole or in part, without written consent of Michael Farrar. For
- further information regarding permission for use or reproduction, please
- contact: Michael Farrar at farrar.michael@gmail.com.
-*******************************************************************/
-
-/*
- Written by Michael Farrar, 2006.
- Please send bug reports and/or suggestions to farrar.michael@gmail.com.
-*/
-
-#include <stdio.h>
-
-#include "defs.h"
-#include "param.h"
-#include "dropgsw.h"
-#include "smith_waterman_sse2.h"
-
-#ifdef __SUNPRO_C
-#include <sunmedia_intrin.h>
-#else
-#include <emmintrin.h>
-#endif
-
-#ifdef SW_SSE2
-
-int
-smith_waterman_sse2_word(const unsigned char * query_sequence,
- unsigned short * query_profile_word,
- const int query_length,
- const unsigned char * db_sequence,
- const int db_length,
- unsigned short gap_open,
- unsigned short gap_extend,
- struct f_struct * f_str)
-{
- int i, j, k;
- short score;
-
- int cmp;
- int iter = (query_length + 7) / 8;
-
- __m128i *p;
- __m128i *workspace = (__m128i *) f_str->workspace;
-
- __m128i E, F, H;
-
- __m128i v_maxscore;
- __m128i v_gapopen;
- __m128i v_gapextend;
-
- __m128i v_min;
- __m128i v_minimums;
- __m128i v_temp;
-
- __m128i *pHLoad, *pHStore;
- __m128i *pE;
-
- __m128i *pScore;
-
- /* Load gap opening penalty to all elements of a constant */
- v_gapopen = _mm_insert_epi16 (v_gapopen, gap_open, 0);
- v_gapopen = _mm_shufflelo_epi16 (v_gapopen, 0);
- v_gapopen = _mm_shuffle_epi32 (v_gapopen, 0);
-
- /* Load gap extension penalty to all elements of a constant */
- v_gapextend = _mm_insert_epi16 (v_gapextend, gap_extend, 0);
- v_gapextend = _mm_shufflelo_epi16 (v_gapextend, 0);
- v_gapextend = _mm_shuffle_epi32 (v_gapextend, 0);
-
- /* load v_maxscore with the zeros. since we are using signed */
- /* math, we will bias the maxscore to -32768 so we have the */
- /* full range of the short. */
- v_maxscore = _mm_cmpeq_epi16 (v_maxscore, v_maxscore);
- v_maxscore = _mm_slli_epi16 (v_maxscore, 15);
-
- v_minimums = _mm_shuffle_epi32 (v_maxscore, 0);
-
- v_min = _mm_shuffle_epi32 (v_maxscore, 0);
- v_min = _mm_srli_si128 (v_min, 14);
-
- /* Zero out the storage vector */
- k = 2 * iter;
-
- p = workspace;
- for (i = 0; i < k; i++)
- {
- _mm_store_si128 (p++, v_maxscore);
- }
-
- pE = workspace;
- pHStore = pE + iter;
- pHLoad = pHStore + iter;
-
- for (i = 0; i < db_length; ++i)
- {
- /* fetch first data asap. */
- pScore = (__m128i *) query_profile_word + db_sequence[i] * iter;
-
- /* bias all elements in F to -32768 */
- F = _mm_cmpeq_epi16 (F, F);
- F = _mm_slli_epi16 (F, 15);
-
- /* load the next h value */
- H = _mm_load_si128 (pHStore + iter - 1);
- H = _mm_slli_si128 (H, 2);
- H = _mm_or_si128 (H, v_min);
-
- p = pHLoad;
- pHLoad = pHStore;
- pHStore = p;
-
- for (j = 0; j < iter; j++)
- {
- /* load E values */
- E = _mm_load_si128 (pE + j);
-
- /* add score to H */
- H = _mm_adds_epi16 (H, *pScore++);
-
- /* Update highest score encountered this far */
- v_maxscore = _mm_max_epi16 (v_maxscore, H);
-
- /* get max from H, E and F */
- H = _mm_max_epi16 (H, E);
- H = _mm_max_epi16 (H, F);
-
- /* save H values */
- _mm_store_si128 (pHStore + j, H);
-
- /* subtract the gap open penalty from H */
- H = _mm_subs_epi16 (H, v_gapopen);
-
- /* update E value */
- E = _mm_subs_epi16 (E, v_gapextend);
- E = _mm_max_epi16 (E, H);
-
- /* update F value */
- F = _mm_subs_epi16 (F, v_gapextend);
- F = _mm_max_epi16 (F, H);
-
- /* save E values */
- _mm_store_si128 (pE + j, E);
-
- /* load the next h value */
- H = _mm_load_si128 (pHLoad + j);
- }
-
- /* reset pointers to the start of the saved data */
- j = 0;
- H = _mm_load_si128 (pHStore + j);
-
- /* the computed F value is for the given column. since */
- /* we are at the end, we need to shift the F value over */
- /* to the next column. */
- F = _mm_slli_si128 (F, 2);
- F = _mm_or_si128 (F, v_min);
- v_temp = _mm_subs_epi16 (H, v_gapopen);
- v_temp = _mm_cmpgt_epi16 (F, v_temp);
- cmp = _mm_movemask_epi8 (v_temp);
-
- while (cmp != 0x0000)
- {
- E = _mm_load_si128 (pE + j);
-
- H = _mm_max_epi16 (H, F);
-
- /* save H values */
- _mm_store_si128 (pHStore + j, H);
-
- /* update E in case the new H value would change it */
- H = _mm_subs_epi16 (H, v_gapopen);
- E = _mm_max_epi16 (E, H);
- _mm_store_si128 (pE + j, E);
-
- /* update F value */
- F = _mm_subs_epi16 (F, v_gapextend);
-
- j++;
- if (j >= iter)
- {
- j = 0;
- F = _mm_slli_si128 (F, 2);
- F = _mm_or_si128 (F, v_min);
- }
- H = _mm_load_si128 (pHStore + j);
-
- v_temp = _mm_subs_epi16 (H, v_gapopen);
- v_temp = _mm_cmpgt_epi16 (F, v_temp);
- cmp = _mm_movemask_epi8 (v_temp);
- }
- }
-
- /* find largest score in the v_maxscore vector */
- v_temp = _mm_srli_si128 (v_maxscore, 8);
- v_maxscore = _mm_max_epi16 (v_maxscore, v_temp);
- v_temp = _mm_srli_si128 (v_maxscore, 4);
- v_maxscore = _mm_max_epi16 (v_maxscore, v_temp);
- v_temp = _mm_srli_si128 (v_maxscore, 2);
- v_maxscore = _mm_max_epi16 (v_maxscore, v_temp);
-
- /* extract the largest score */
- score = _mm_extract_epi16 (v_maxscore, 0);
-
- /* return largest score biased by 32768 */
- return score + 32768;
-}
-
-
-
-
-int
-smith_waterman_sse2_byte(const unsigned char * query_sequence,
- unsigned char * query_profile_byte,
- const int query_length,
- const unsigned char * db_sequence,
- const int db_length,
- unsigned char bias,
- unsigned char gap_open,
- unsigned char gap_extend,
- struct f_struct * f_str)
-{
- int i, j, k;
- int score;
-
- int dup;
- int cmp;
- int iter = (query_length + 15) / 16;
-
- __m128i *p;
- __m128i *workspace = (__m128i *) f_str->workspace;
-
- __m128i E, F, H;
-
- __m128i v_maxscore;
- __m128i v_bias;
- __m128i v_gapopen;
- __m128i v_gapextend;
-
- __m128i v_temp;
- __m128i v_zero;
-
- __m128i *pHLoad, *pHStore;
- __m128i *pE;
-
- __m128i *pScore;
-
- /* Load the bias to all elements of a constant */
- dup = ((short) bias << 8) | bias;
- v_bias = _mm_insert_epi16 (v_bias, dup, 0);
- v_bias = _mm_shufflelo_epi16 (v_bias, 0);
- v_bias = _mm_shuffle_epi32 (v_bias, 0);
-
- /* Load gap opening penalty to all elements of a constant */
- dup = ((short) gap_open << 8) | gap_open;
- v_gapopen = _mm_insert_epi16 (v_gapopen, dup, 0);
- v_gapopen = _mm_shufflelo_epi16 (v_gapopen, 0);
- v_gapopen = _mm_shuffle_epi32 (v_gapopen, 0);
-
- /* Load gap extension penalty to all elements of a constant */
- dup = ((short) gap_extend << 8) | gap_extend;
- v_gapextend = _mm_insert_epi16 (v_gapextend, dup, 0);
- v_gapextend = _mm_shufflelo_epi16 (v_gapextend, 0);
- v_gapextend = _mm_shuffle_epi32 (v_gapextend, 0);
-
- /* initialize the max score */
- v_maxscore = _mm_xor_si128 (v_maxscore, v_maxscore);
-
- /* create a constant of all zeros for comparison */
- v_zero = _mm_xor_si128 (v_zero, v_zero);
-
- /* Zero out the storage vector */
- k = iter * 2;
-
- p = workspace;
- for (i = 0; i < k; i++)
- {
- _mm_store_si128 (p++, v_maxscore);
- }
-
- pE = workspace;
- pHStore = pE + iter;
- pHLoad = pHStore + iter;
-
- for (i = 0; i < db_length; ++i)
- {
- /* fetch first data asap. */
- pScore = (__m128i *) query_profile_byte + db_sequence[i] * iter;
-
- /* zero out F value. */
- F = _mm_xor_si128 (F, F);
-
- /* load the next h value */
- H = _mm_load_si128 (pHStore + iter - 1);
- H = _mm_slli_si128 (H, 1);
-
- p = pHLoad;
- pHLoad = pHStore;
- pHStore = p;
-
- for (j = 0; j < iter; j++)
- {
- /* load values E. */
- E = _mm_load_si128 (pE + j);
-
- /* add score to H */
- H = _mm_adds_epu8 (H, *pScore++);
- H = _mm_subs_epu8 (H, v_bias);
-
- /* Update highest score encountered this far */
- v_maxscore = _mm_max_epu8 (v_maxscore, H);
-
- /* get max from H, E and F */
- H = _mm_max_epu8 (H, E);
- H = _mm_max_epu8 (H, F);
-
- /* save H values */
- _mm_store_si128 (pHStore + j, H);
-
- /* subtract the gap open penalty from H */
- H = _mm_subs_epu8 (H, v_gapopen);
-
- /* update E value */
- E = _mm_subs_epu8 (E, v_gapextend);
- E = _mm_max_epu8 (E, H);
-
- /* update F value */
- F = _mm_subs_epu8 (F, v_gapextend);
- F = _mm_max_epu8 (F, H);
-
- /* save E values */
- _mm_store_si128 (pE + j, E);
-
- /* load the next h value */
- H = _mm_load_si128 (pHLoad + j);
- }
-
- /* reset pointers to the start of the saved data */
- j = 0;
- H = _mm_load_si128 (pHStore + j);
-
- /* the computed F value is for the given column. since */
- /* we are at the end, we need to shift the F value over */
- /* to the next column. */
- F = _mm_slli_si128 (F, 1);
- v_temp = _mm_subs_epu8 (H, v_gapopen);
- v_temp = _mm_subs_epu8 (F, v_temp);
- v_temp = _mm_cmpeq_epi8 (v_temp, v_zero);
- cmp = _mm_movemask_epi8 (v_temp);
-
- while (cmp != 0xffff)
- {
- E = _mm_load_si128 (pE + j);
-
- H = _mm_max_epu8 (H, F);
-
- /* save H values */
- _mm_store_si128 (pHStore + j, H);
-
- /* update E in case the new H value would change it */
- H = _mm_subs_epu8 (H, v_gapopen);
- E = _mm_max_epu8 (E, H);
- _mm_store_si128 (pE + j, E);
-
- /* update F value */
- F = _mm_subs_epu8 (F, v_gapextend);
-
- j++;
- if (j >= iter)
- {
- j = 0;
- F = _mm_slli_si128 (F, 1);
- }
- H = _mm_load_si128 (pHStore + j);
-
- v_temp = _mm_subs_epu8 (H, v_gapopen);
- v_temp = _mm_subs_epu8 (F, v_temp);
- v_temp = _mm_cmpeq_epi8 (v_temp, v_zero);
- cmp = _mm_movemask_epi8 (v_temp);
- }
- }
-
- /* find largest score in the v_maxscore vector */
- v_temp = _mm_srli_si128 (v_maxscore, 8);
- v_maxscore = _mm_max_epu8 (v_maxscore, v_temp);
- v_temp = _mm_srli_si128 (v_maxscore, 4);
- v_maxscore = _mm_max_epu8 (v_maxscore, v_temp);
- v_temp = _mm_srli_si128 (v_maxscore, 2);
- v_maxscore = _mm_max_epu8 (v_maxscore, v_temp);
- v_temp = _mm_srli_si128 (v_maxscore, 1);
- v_maxscore = _mm_max_epu8 (v_maxscore, v_temp);
-
- /* store in temporary variable */
- score = _mm_extract_epi16 (v_maxscore, 0);
- score = score & 0x00ff;
-
- /* check if we might have overflowed */
- if (score + bias >= 255)
- {
- score = 255;
- }
-
- /* return largest score */
- return score;
-}
-#else
-
-/* No SSE2 support. Avoid compiler complaints about empty object */
-
-int sw_dummy;
-
-#endif
git://source.jalview.org / jabaws.git / blobdiff