diagonal_matrix.h

Go to the documentation of this file.

template<class T, class S> class diagonal_matrix : public matrix<T, S>
{
private:
        vector<S> _val;         
public:
        diagonal_matrix(const vector<S> &_val)
        {
                this->_val = _val;
        }
        void multiply(const packed_vector<S> &_x, packed_vector<S> &_b) const
        {
                const S *a = &_val[0];
                const S *x = &_x[0];
                S *b = &_b[0];
                int m = (int)_val.size();
                int n = _x.numrhs();
                switch(n)
                {
                case 0:
                        break;
                case 1:
                        for(int i = 0; i < m; i++)
                        {
                                S val = a[i];
                                b[i] = val * x[i];
                        }
                        break;
                case 2:
                        for(int i = 0; i < m; i++)
                        {
                                T k = 2 * i;
                                S val = a[i];
                                b[k + 0] = val * x[k + 0];
                                b[k + 1] = val * x[k + 1];
                        }
                        break;
                case 3:
                        for(int i = 0; i < m; i++)
                        {
                                T k = 3 * i;
                                S val = a[i];
                                b[k + 0] = val * x[k + 0];
                                b[k + 1] = val * x[k + 1];
                                b[k + 2] = val * x[k + 2];
                        }
                        break;
                case 4:
                        for(int i = 0; i < m; i++)
                        {
                                T k = 4 * i;
                                S val = a[i];
                                b[k + 0] = val * x[k + 0];
                                b[k + 1] = val * x[k + 1];
                                b[k + 2] = val * x[k + 2];
                                b[k + 3] = val * x[k + 3];
                        }
                        break;
                default:
                        for(int i = 0; i < m; i++)
                        {
                                T k = n * i;
                                S val = a[i];
                                for(int j = 0; j < n; j++)
                                {
                                        b[k + j] = val * x[k + j];
                                }
                        }
                }
        }
};

---