1 /*
2 Temporary implementation of vector and matrix classes
3 No attempt is made to check if the function arguments are valid
4 */
5
6 #include <iostream>
7 #include <math.h> // for sqrt()
8
9 #include "linalg.h"
10
eps()11 double eps() {
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 ;
15 Current = 1.0 ;
16 do
17 {
18 Last = Current ;
19 Current /= 2.0 ;
20 OnePlusCurrent = 1.0 + Current ;
21 } while (OnePlusCurrent > 1.0) ;
22 return Last ;
23 }
24
25
RVector()26 RVector::RVector() {
27 // Constructor
28 len=0;
29 elements=0; (*this)=0.;
30 }
31
RVector(int n)32 RVector::RVector(int n) {
33 // Constructor
34 len=n;
35 elements=new double[len]; (*this)=0.;
36 }
37
RVector(RCRVector vect)38 RVector::RVector(RCRVector vect)
39 {
40 // Constructor + Copy
41 len=vect.len;
42 elements=new double[len]; (*this)=vect;
43 }
44
operator =(RCRVector vect)45 RCRVector RVector::operator=(RCRVector vect)
46 {
47 // Copy
48 for (int i=0;i<len;i++) elements[i]=vect.elements[i];
49 return *this;
50 }
51
operator =(double num)52 RCRVector RVector::operator=(double num) {
53 // Assignment
54 for (int i=0;i<len;i++) elements[i]=num;
55 return *this;
56 }
57
nrm2()58 double RVector::nrm2() {
59 double sum=0 ;
60 for (int i = 0; i < len; i++)
61 sum+=elements[i]*elements[i] ;
62 return sqrt(sum) ;
63 }
64
scal(double alpha,RCRVector x)65 void scal(double alpha, RCRVector x) {
66 int n=x.len ;
67 for (int i = 0; i < n; i++)
68 x.elements[i] = alpha*x.elements[i] ;
69 }
70
norm2(RCRVector x)71 double norm2(RCRVector x) {
72 // Euclidian norm
73 double sum=0 ;
74 double* pa=x.elements ;
75 int n=x.len ;
76 for (int i = 0; i < n; i++) {
77 sum+=(*pa)*(*pa) ;
78 pa++ ;
79 }
80 return sqrt(sum) ;
81 }
82
normInf(RCRVector x)83 double normInf(RCRVector x) {
84 // Infinity norm
85 int n=x.len ;
86 double tmp=DBL_MIN ;
87 for (int i=0 ; i<n ; i++)
88 tmp=max(tmp,fabs(x.elements[i])) ;
89 return tmp ;
90 }
91
dot(RCRVector x,RCRVector y)92 double dot(RCRVector x, RCRVector y) {
93 // dot <- x'y
94 int n=x.len ;
95 double sum=0 ;
96 for (int i=0 ; i<n ; i++)
97 sum+= x.elements[i]*y.elements[i] ;
98 return sum ;
99 }
100
copy(RCRVector x,RCRVector y)101 void copy(RCRVector x, RCRVector y) {
102 // y <- x
103 double *px=x.elements ;
104 double *py=y.elements ;
105 int n=x.len ;
106 for (int i=0 ; i<n ; i++)
107 (*py++)=(*px++) ;
108 }
109
axpy(double alpha,RCRVector x,RCRVector y)110 void axpy(double alpha, RCRVector x, RCRVector y) {
111 // y <- alpha*x + y
112 double *px=x.elements ;
113 double *py=y.elements ;
114 int n=x.len ;
115 for (int i=0 ; i<n ; i++) {
116 *py=alpha*(*px) + *py ;
117 px++ ; py++ ;
118 }
119 }
120
gemv(char trans,double alpha,RCRMatrix A,RCRVector x,double beta,RCRVector y)121 void gemv(char trans, double alpha, RCRMatrix A,RCRVector x,
122 double beta, RCRVector y) {
123 // Matrix-vector multiplication
124 int i, j, dim=A.Dim;
125 double sum ;
126
127 if (trans=='N') {
128 // y=alpha*A*x + beta*y
129 for (i=0;i<dim;i++) {
130 sum=0.;
131 for (j=0;j<dim;j++)
132 sum+=A.Vals[j+i*dim]*x.elements[j]*alpha;
133 y.elements[i]=y.elements[i]*beta + sum ;
134 }
135 }
136 else {
137 // y=alpha*transpose(A)*x + beta*y
138 for (i=0;i<dim;i++) {
139 sum=0.;
140 for (j=0;j<dim;j++) {
141 sum+=A.Vals[i+j*dim]*x.elements[j]*alpha ;
142 }
143 y.elements[i]=y.elements[i]*beta + sum ;
144 }
145 }
146 }
147
operator <<(ostream & os,const RVector & v)148 ostream & operator << (ostream & os, const RVector & v) {
149 os << '[';
150 for (int i = 0; i < v.len; i++) {
151 if (i>0) os << "," ;
152 os << v.elements[i] ;
153 }
154 return os << ']';
155 }
156
157 /*************************** Matrix Class ***********************/
158
RMatrix()159 RMatrix::RMatrix() {
160 // Constructor
161 Dim=0 ; Vals=0 ; (*this)=0 ;
162 }
163
RMatrix(int dim)164 RMatrix::RMatrix(int dim) {
165 // Constructor
166 Dim=dim;
167 Vals=new double[long(Dim)*long(Dim)]; (*this)=0.;
168 }
169
RMatrix(RCRMatrix matr)170 RMatrix::RMatrix(RCRMatrix matr) {
171 // Constructor + Copy
172 Dim=matr.Dim;
173 Vals=new double[long(Dim)*long(Dim)]; (*this)=matr;
174 }
175
operator =(RCRMatrix mat)176 RCRMatrix RMatrix::operator=(RCRMatrix mat)
177 { // Assignment, A=B
178 long int Len=long(Dim)*long(Dim);
179 for (int i=0;i<Len;i++) Vals[i]=mat.Vals[i] ;
180 return *this;
181 }
182
operator =(double num)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;
186 return *this;
187 }
188
operator ()(int vidx,int hidx)189 double& RMatrix::operator()(int vidx,int hidx) {
190 return Vals[vidx*Dim+hidx];
191 }
192
ger(double alpha,RCRVector x,RCRVector y,RCRMatrix A)193 void ger(double alpha, RCRVector x,RCRVector y, RCRMatrix A) {
194 // Rank one update : A=alpha*xy'+A
195 int dim=x.len;
196 double* pa=A.Vals ;
197
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;
201 pa++ ;
202 }
203 }
204
operator <<(ostream & os,const RMatrix & A)205 ostream & operator << (ostream & os, const RMatrix & A) {
206 int n=A.Dim ;
207 double* pa=A.Vals ;
208 os << endl ;
209 for (int i = 0; i < n; i++) {
210 for (int j=0 ; j< n ; j++) {
211 os << (*(pa++)) << " " ;
212 }
213 os << endl ;
214 }
215 return os ;
216 }
217