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