stogo/linalg.cc

   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
  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
  26 RVector::RVector() {
  27  // Constructor
  28  len=0;
  29  elements=0; (*this)=0.;
  30 }
  31
  32 RVector::RVector(int n) {
  33  // Constructor
  34  len=n;
  35  elements=new double[len]; (*this)=0.;
  36 }
  37
  38 RVector::RVector(RCRVector vect)
  39 {
  40  // Constructor + Copy
  41  len=vect.len;
  42  elements=new double[len]; (*this)=vect;
  43 }
  44
  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
  52 RCRVector RVector::operator=(double num) {
  53  // Assignment
  54  for (int i=0;i<len;i++) elements[i]=num;
  55  return *this;
  56 }
  57
  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
  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
  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
  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
  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
 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
 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
 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
 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
 159 RMatrix::RMatrix() {
 160  // Constructor
 161  Dim=0 ; Vals=0 ; (*this)=0 ;
 162 }
 163
 164 RMatrix::RMatrix(int dim) {
 165  // Constructor
 166  Dim=dim;
 167  Vals=new double[long(Dim)*long(Dim)]; (*this)=0.;
 168 }
 169
 170 RMatrix::RMatrix(RCRMatrix matr) {
 171  // Constructor + Copy
 172  Dim=matr.Dim;
 173  Vals=new double[long(Dim)*long(Dim)]; (*this)=matr;
 174 }
 175
 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
 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
 189 double& RMatrix::operator()(int vidx,int hidx) {
 190  return Vals[vidx*Dim+hidx];
 191 }
 192
 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
 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 }