2 Temporary implementation of vector and matrix classes
3 No attempt is made to check if the function arguments are valid
7 #include <math.h> // for sqrt()
12 /* Returns the machine precision : (min { x >= 0 : 1 + x > 1 })
13 NB This routine should be replaced by LAPACK_LAMCH */
14 double Current, Last, OnePlusCurrent ;
20 OnePlusCurrent = 1.0 + Current ;
21 } while (OnePlusCurrent > 1.0) ;
29 elements=0; (*this)=0.;
32 RVector::RVector(int n) {
35 elements=new double[len]; (*this)=0.;
38 RVector::RVector(RCRVector vect)
42 elements=new double[len]; (*this)=vect;
45 RCRVector RVector::operator=(RCRVector vect)
48 for (int i=0;i<len;i++) elements[i]=vect.elements[i];
52 RCRVector RVector::operator=(double num) {
54 for (int i=0;i<len;i++) elements[i]=num;
58 double RVector::nrm2() {
60 for (int i = 0; i < len; i++)
61 sum+=elements[i]*elements[i] ;
65 void scal(double alpha, RCRVector x) {
67 for (int i = 0; i < n; i++)
68 x.elements[i] = alpha*x.elements[i] ;
71 double norm2(RCRVector x) {
74 double* pa=x.elements ;
76 for (int i = 0; i < n; i++) {
83 double normInf(RCRVector x) {
87 for (int i=0 ; i<n ; i++)
88 tmp=max(tmp,fabs(x.elements[i])) ;
92 double dot(RCRVector x, RCRVector y) {
96 for (int i=0 ; i<n ; i++)
97 sum+= x.elements[i]*y.elements[i] ;
101 void copy(RCRVector x, RCRVector y) {
103 double *px=x.elements ;
104 double *py=y.elements ;
106 for (int i=0 ; i<n ; i++)
110 void axpy(double alpha, RCRVector x, RCRVector y) {
112 double *px=x.elements ;
113 double *py=y.elements ;
115 for (int i=0 ; i<n ; i++) {
116 *py=alpha*(*px) + *py ;
121 void gemv(char trans, double alpha, RCRMatrix A,RCRVector x,
122 double beta, RCRVector y) {
123 // Matrix-vector multiplication
128 // y=alpha*A*x + beta*y
129 for (i=0;i<dim;i++) {
132 sum+=A.Vals[j+i*dim]*x.elements[j]*alpha;
133 y.elements[i]=y.elements[i]*beta + sum ;
137 // y=alpha*transpose(A)*x + beta*y
138 for (i=0;i<dim;i++) {
140 for (j=0;j<dim;j++) {
141 sum+=A.Vals[i+j*dim]*x.elements[j]*alpha ;
143 y.elements[i]=y.elements[i]*beta + sum ;
148 ostream & operator << (ostream & os, const RVector & v) {
150 for (int i = 0; i < v.len; i++) {
152 os << v.elements[i] ;
157 /*************************** Matrix Class ***********************/
161 Dim=0 ; Vals=0 ; (*this)=0 ;
164 RMatrix::RMatrix(int dim) {
167 Vals=new double[long(Dim)*long(Dim)]; (*this)=0.;
170 RMatrix::RMatrix(RCRMatrix matr) {
171 // Constructor + Copy
173 Vals=new double[long(Dim)*long(Dim)]; (*this)=matr;
176 RCRMatrix RMatrix::operator=(RCRMatrix mat)
178 long int Len=long(Dim)*long(Dim);
179 for (int i=0;i<Len;i++) Vals[i]=mat.Vals[i] ;
183 RCRMatrix RMatrix::operator=(double num) {
184 long int Len=long(Dim)*long(Dim);
185 for (long int i=0;i<Len;i++) Vals[i]=num;
189 double& RMatrix::operator()(int vidx,int hidx) {
190 return Vals[vidx*Dim+hidx];
193 void ger(double alpha, RCRVector x,RCRVector y, RCRMatrix A) {
194 // Rank one update : A=alpha*xy'+A
198 for (int i=0;i<dim;i++)
199 for (int j=0;j<dim;j++) {
200 *pa=alpha*x.elements[i]*y.elements[j] + *pa;
205 ostream & operator << (ostream & os, const RMatrix & A) {
209 for (int i = 0; i < n; i++) {
210 for (int j=0 ; j< n ; j++) {
211 os << (*(pa++)) << " " ;