| 28 | #include "sblas.h" |
| 29 | |
| 30 | void symmetricmatrixvectormultiply(const ap::real_2d_array& a, |
| 31 | bool isupper, |
| 32 | int i1, |
| 33 | int i2, |
| 34 | const ap::real_1d_array& x, |
| 35 | double alpha, |
| 36 | ap::real_1d_array& y) |
| 37 | { |
| 38 | int i; |
| 39 | int ba1; |
| 40 | int ba2; |
| 41 | int by1; |
| 42 | int by2; |
| 43 | int bx1; |
| 44 | int bx2; |
| 45 | int n; |
| 46 | double v; |
| 47 | |
| 48 | n = i2-i1+1; |
| 49 | if( n<=0 ) |
| 50 | { |
| 51 | return; |
| 52 | } |
| 53 | |
| 54 | // |
| 55 | // Let A = L + D + U, where |
| 56 | // L is strictly lower triangular (main diagonal is zero) |
| 57 | // D is diagonal |
| 58 | // U is strictly upper triangular (main diagonal is zero) |
| 59 | // |
| 60 | // A*x = L*x + D*x + U*x |
| 61 | // |
| 62 | // Calculate D*x first |
| 63 | // |
| 64 | for(i = i1; i <= i2; i++) |
| 65 | { |
| 66 | y(i-i1+1) = a(i,i)*x(i-i1+1); |
| 67 | } |
| 68 | |
| 69 | // |
| 70 | // Add L*x + U*x |
| 71 | // |
| 72 | if( isupper ) |
| 73 | { |
| 74 | for(i = i1; i <= i2-1; i++) |
| 75 | { |
| 76 | |
| 77 | // |
| 78 | // Add L*x to the result |
| 79 | // |
| 80 | v = x(i-i1+1); |
| 81 | by1 = i-i1+2; |
| 82 | by2 = n; |
| 83 | ba1 = i+1; |
| 84 | ba2 = i2; |
| 85 | ap::vadd(&y(by1), 1, &a(i, ba1), 1, ap::vlen(by1,by2), v); |
| 86 | |
| 87 | // |
no test coverage detected