00001
00005 template<class T, class S> class idx_matrix : public matrix<T, S>
00006 {
00007 private:
00008 vector<T> _row;
00009 vector<T> _col;
00010 vector<S> _val;
00011 public:
00018 idx_matrix(const vector<T> &_row, const vector<T> &_col, const vector<S> &_val)
00019 {
00020 this->_row = _row;
00021 this->_col = _col;
00022 this->_val = _val;
00023 }
00029 void multiply(const packed_vector<S> &_x, packed_vector<S> &_b) const
00030 {
00031 memset(&_b[0], 0, sizeof(S) * _b.size());
00032 const T *row = &_row[0], *col = &_col[0];
00033 const S *a = &_val[0];
00034 const S *x = &_x[0];
00035 S *b = &_b[0];
00036 int m = (int)_val.size();
00037 int n = _x.numrhs();
00038 switch(n)
00039 {
00040 case 0:
00041 break;
00042 case 1:
00043 for(int i = 0; i < m; i++)
00044 {
00045 T k = row[i];
00046 T l = col[i];
00047 S val = a[i];
00048 b[k] += val * x[l];
00049 }
00050 break;
00051 case 2:
00052 for(int i = 0; i < m; i++)
00053 {
00054 T k = 2 * row[i];
00055 T l = 2 * col[i];
00056 S val = a[i];
00057 b[k + 0] += val * x[l + 0];
00058 b[k + 1] += val * x[l + 1];
00059 }
00060 break;
00061 case 3:
00062 for(int i = 0; i < m; i++)
00063 {
00064 T k = 3 * row[i];
00065 T l = 3 * col[i];
00066 S val = a[i];
00067 b[k + 0] += val * x[l + 0];
00068 b[k + 1] += val * x[l + 1];
00069 b[k + 2] += val * x[l + 2];
00070 }
00071 break;
00072 case 4:
00073 for(int i = 0; i < m; i++)
00074 {
00075 T k = 4 * row[i];
00076 T l = 4 * col[i];
00077 S val = a[i];
00078 b[k + 0] += val * x[l + 0];
00079 b[k + 1] += val * x[l + 1];
00080 b[k + 2] += val * x[l + 2];
00081 b[k + 3] += val * x[l + 3];
00082 }
00083 break;
00084 case 8:
00085 for(int i = 0; i < m; i++)
00086 {
00087 T k = 8 * row[i];
00088 T l = 8 * col[i];
00089 S val = a[i];
00090 b[k + 0] += val * x[l + 0];
00091 b[k + 1] += val * x[l + 1];
00092 b[k + 2] += val * x[l + 2];
00093 b[k + 3] += val * x[l + 3];
00094 b[k + 4] += val * x[l + 4];
00095 b[k + 5] += val * x[l + 5];
00096 b[k + 6] += val * x[l + 6];
00097 b[k + 7] += val * x[l + 7];
00098 }
00099 break;
00100 case 16:
00101 for(int i = 0; i < m; i++)
00102 {
00103 T k = 16 * row[i];
00104 T l = 16 * col[i];
00105 S val = a[i];
00106 b[k + 0] += val * x[l + 0];
00107 b[k + 1] += val * x[l + 1];
00108 b[k + 2] += val * x[l + 2];
00109 b[k + 3] += val * x[l + 3];
00110 b[k + 4] += val * x[l + 4];
00111 b[k + 5] += val * x[l + 5];
00112 b[k + 6] += val * x[l + 6];
00113 b[k + 7] += val * x[l + 7];
00114 b[k + 8] += val * x[l + 8];
00115 b[k + 9] += val * x[l + 9];
00116 b[k + 10] += val * x[l + 10];
00117 b[k + 11] += val * x[l + 11];
00118 b[k + 12] += val * x[l + 12];
00119 b[k + 13] += val * x[l + 13];
00120 b[k + 14] += val * x[l + 14];
00121 b[k + 15] += val * x[l + 15];
00122 }
00123 break;
00124 default:
00125 for(int i = 0; i < m; i++)
00126 {
00127 T k = n * row[i];
00128 T l = n * col[i];
00129 S val = a[i];
00130 for(int j = 0; j < n; j++)
00131 {
00132 b[k + j] += val * x[l + j];
00133 }
00134 }
00135 }
00136 }
00137 };
