SNAP Library 4.1, User Reference  2018-07-26 16:30:42
SNAP, a general purpose, high performance system for analysis and manipulation of large networks
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
TNumericalStuff Class Reference

#include <linalg.h>

Static Public Member Functions

static void SymetricToTridiag (TFltVV &a, int n, TFltV &d, TFltV &e)
 
static void EigSymmetricTridiag (TFltV &d, TFltV &e, int n, TFltVV &z)
 
static void CholeskyDecomposition (TFltVV &A, TFltV &p)
 
static void CholeskySolve (const TFltVV &A, const TFltV &p, const TFltV &b, TFltV &x)
 
static void SolveSymetricSystem (TFltVV &A, const TFltV &b, TFltV &x)
 
static void InverseSubstitute (TFltVV &A, const TFltV &p)
 
static void InverseSymetric (TFltVV &A)
 
static void InverseTriagonal (TFltVV &A)
 
static void LUDecomposition (TFltVV &A, TIntV &indx, double &d)
 
static void LUSolve (const TFltVV &A, const TIntV &indx, TFltV &b)
 
static void SolveLinearSystem (TFltVV &A, const TFltV &b, TFltV &x)
 

Static Private Member Functions

static double sqr (double a)
 
static double sign (double a, double b)
 
static double pythag (double a, double b)
 
static void nrerror (const TStr &error_text)
 

Detailed Description

Definition at line 300 of file linalg.h.

Member Function Documentation

void TNumericalStuff::CholeskyDecomposition ( TFltVV A,
TFltV p 
)
static

Definition at line 797 of file linalg.cpp.

797  {
798  Assert(A.GetRows() == A.GetCols());
799  int n = A.GetRows(); p.Reserve(n,n);
800 
801  int i,j,k;
802  double sum;
803  for (i=1;i<=n;i++) {
804  for (j=i;j<=n;j++) {
805  for (sum=A(i-1,j-1),k=i-1;k>=1;k--) sum -= A(i-1,k-1)*A(j-1,k-1);
806  if (i == j) {
807  if (sum <= 0.0)
808  nrerror("choldc failed");
809  p[i-1]=sqrt(sum);
810  } else A(j-1,i-1)=sum/p[i-1];
811  }
812  }
813 }
#define Assert(Cond)
Definition: bd.h:251
static void nrerror(const TStr &error_text)
Definition: linalg.cpp:655
TSizeTy GetRows() const
Definition: ds.h:2252
void Reserve(const TSizeTy &_MxVals)
Reserves enough memory for the vector to store _MxVals elements.
Definition: ds.h:543
TSizeTy GetCols() const
Definition: ds.h:2253
void TNumericalStuff::CholeskySolve ( const TFltVV A,
const TFltV p,
const TFltV b,
TFltV x 
)
static

Definition at line 815 of file linalg.cpp.

815  {
816  IAssert(A.GetRows() == A.GetCols());
817  int n = A.GetRows(); x.Reserve(n,n);
818 
819  int i,k;
820  double sum;
821 
822  // Solve L * y = b, storing y in x
823  for (i=1;i<=n;i++) {
824  for (sum=b[i-1],k=i-1;k>=1;k--)
825  sum -= A(i-1,k-1)*x[k-1];
826  x[i-1]=sum/p[i-1];
827  }
828 
829  // Solve L^T * x = y
830  for (i=n;i>=1;i--) {
831  for (sum=x[i-1],k=i+1;k<=n;k++)
832  sum -= A(k-1,i-1)*x[k-1];
833  x[i-1]=sum/p[i-1];
834  }
835 }
#define IAssert(Cond)
Definition: bd.h:262
TSizeTy GetRows() const
Definition: ds.h:2252
void Reserve(const TSizeTy &_MxVals)
Reserves enough memory for the vector to store _MxVals elements.
Definition: ds.h:543
TSizeTy GetCols() const
Definition: ds.h:2253
void TNumericalStuff::EigSymmetricTridiag ( TFltV d,
TFltV e,
int  n,
TFltVV z 
)
static

Definition at line 740 of file linalg.cpp.

740  {
741  int m,l,iter,i,k; // N = n+1;
742  double s,r,p,g,f,dd,c,b;
743  // Convenient to renumber the elements of e
744  for (i=2;i<=n;i++) e[i-1]=e[i];
745  e[n]=0.0;
746  for (l=1;l<=n;l++) {
747  iter=0;
748  do {
749  // Look for a single small subdiagonal element to split the matrix.
750  for (m=l;m<=n-1;m++) {
751  dd=TFlt::Abs(d[m])+TFlt::Abs(d[m+1]);
752  if ((double)(TFlt::Abs(e[m])+dd) == dd) break;
753  }
754  if (m != l) {
755  if (iter++ == 60) nrerror("Too many iterations in EigSymmetricTridiag");
756  //Form shift.
757  g=(d[l+1]-d[l])/(2.0*e[l]);
758  r=pythag(g,1.0);
759  //This is dm - ks.
760  g=d[m]-d[l]+e[l]/(g+sign(r,g));
761  s=c=1.0;
762  p=0.0;
763  // A plane rotation as in the original QL, followed by
764  // Givens rotations to restore tridiagonal form
765  for (i=m-1;i>=l;i--) {
766  f=s*e[i];
767  b=c*e[i];
768  e[i+1]=(r=pythag(f,g));
769  // Recover from underflow.
770  if (r == 0.0) {
771  d[i+1] -= p;
772  e[m]=0.0;
773  break;
774  }
775  s=f/r;
776  c=g/r;
777  g=d[i+1]-p;
778  r=(d[i]-g)*s+2.0*c*b;
779  d[i+1]=g+(p=s*r);
780  g=c*r-b;
781  // Next loop can be omitted if eigenvectors not wanted
782  for (k=0;k<n;k++) {
783  f=z(k,i);
784  z(k,i)=s*z(k,i-1)+c*f;
785  z(k,i-1)=c*z(k,i-1)-s*f;
786  }
787  }
788  if (r == 0.0 && i >= l) continue;
789  d[l] -= p;
790  e[l]=g;
791  e[m]=0.0;
792  }
793  } while (m != l);
794  }
795 }
static void nrerror(const TStr &error_text)
Definition: linalg.cpp:655
static double Abs(const double &Flt)
Definition: dt.h:1427
static double pythag(double a, double b)
Definition: linalg.cpp:660
static double sign(double a, double b)
Definition: linalg.cpp:651
void TNumericalStuff::InverseSubstitute ( TFltVV A,
const TFltV p 
)
static

Definition at line 843 of file linalg.cpp.

843  {
844  IAssert(A.GetRows() == A.GetCols());
845  int n = A.GetRows(); TFltV x(n);
846 
847  int i, j, k; double sum;
848  for (i = 0; i < n; i++) {
849  // solve L * y = e_i, store in x
850  // elements from 0 to i-1 are 0.0
851  for (j = 0; j < i; j++) x[j] = 0.0;
852  // solve l_ii * y_i = 1 => y_i = 1/l_ii
853  x[i] = 1/p[i];
854  // solve y_j for j > i
855  for (j = i+1; j < n; j++) {
856  for (sum = 0.0, k = i; k < j; k++)
857  sum -= A(j,k) * x[k];
858  x[j] = sum / p[j];
859  }
860 
861  // solve L'* x = y, store in upper triangule of A
862  for (j = n-1; j >= i; j--) {
863  for (sum = x[j], k = j+1; k < n; k++)
864  sum -= A(k,j)*x[k];
865  x[j] = sum/p[j];
866  }
867  for (int j = i; j < n; j++) A(i,j) = x[j];
868  }
869 
870 }
#define IAssert(Cond)
Definition: bd.h:262
TSizeTy GetRows() const
Definition: ds.h:2252
TSizeTy GetCols() const
Definition: ds.h:2253
void TNumericalStuff::InverseSymetric ( TFltVV A)
static

Definition at line 872 of file linalg.cpp.

872  {
873  IAssert(A.GetRows() == A.GetCols());
874  TFltV p;
875  // first we calculate cholesky decomposition of A
876  CholeskyDecomposition(A, p);
877  // than we solve system A x_i = e_i for i = 1..n
878  InverseSubstitute(A, p);
879 }
#define IAssert(Cond)
Definition: bd.h:262
static void CholeskyDecomposition(TFltVV &A, TFltV &p)
Definition: linalg.cpp:797
TSizeTy GetRows() const
Definition: ds.h:2252
static void InverseSubstitute(TFltVV &A, const TFltV &p)
Definition: linalg.cpp:843
TSizeTy GetCols() const
Definition: ds.h:2253
void TNumericalStuff::InverseTriagonal ( TFltVV A)
static

Definition at line 881 of file linalg.cpp.

881  {
882  IAssert(A.GetRows() == A.GetCols());
883  int n = A.GetRows(); TFltV x(n), p(n);
884 
885  int i, j, k; double sum;
886  // copy upper triangle to lower one as we'll overwrite upper one
887  for (i = 0; i < n; i++) {
888  p[i] = A(i,i);
889  for (j = i+1; j < n; j++)
890  A(j,i) = A(i,j);
891  }
892  // solve
893  for (i = 0; i < n; i++) {
894  // solve R * x = e_i, store in x
895  // elements from 0 to i-1 are 0.0
896  for (j = n-1; j > i; j--) x[j] = 0.0;
897  // solve l_ii * y_i = 1 => y_i = 1/l_ii
898  x[i] = 1/p[i];
899  // solve y_j for j > i
900  for (j = i-1; j >= 0; j--) {
901  for (sum = 0.0, k = i; k > j; k--)
902  sum -= A(k,j) * x[k];
903  x[j] = sum / p[j];
904  }
905  for (int j = 0; j <= i; j++) A(j,i) = x[j];
906  }
907 }
#define IAssert(Cond)
Definition: bd.h:262
TSizeTy GetRows() const
Definition: ds.h:2252
TSizeTy GetCols() const
Definition: ds.h:2253
void TNumericalStuff::LUDecomposition ( TFltVV A,
TIntV indx,
double &  d 
)
static

Definition at line 909 of file linalg.cpp.

909  {
910  Assert(A.GetRows() == A.GetCols());
911  int n = A.GetRows(); indx.Reserve(n,n);
912 
913  int i=0,imax=0,j=0,k=0;
914  double big,dum,sum,temp;
915  TFltV vv(n); // vv stores the implicit scaling of each row.
916  d=1.0; // No row interchanges yet.
917 
918  // Loop over rows to get the implicit scaling information.
919  for (i=1;i<=n;i++) {
920  big=0.0;
921  for (j=1;j<=n;j++)
922  if ((temp=TFlt::Abs(A(i-1,j-1))) > big) big=temp;
923  if (big == 0.0) nrerror("Singular matrix in routine LUDecomposition");
924  vv[i-1]=1.0/big;
925  }
926 
927  for (j=1;j<=n;j++) {
928  for (i=1;i<j;i++) {
929  sum=A(i-1,j-1);
930  for (k=1;k<i;k++) sum -= A(i-1,k-1)*A(k-1,j-1);
931  A(i-1,j-1)=sum;
932  }
933  big=0.0; //Initialize for the search for largest pivot element.
934  for (i=j;i<=n;i++) {
935  sum=A(i-1,j-1);
936  for (k=1;k<j;k++)
937  sum -= A(i-1,k-1)*A(k-1,j-1);
938  A(i-1,j-1)=sum;
939 
940  //Is the figure of merit for the pivot better than the best so far?
941  if ((dum=vv[i-1] * TFlt::Abs(sum)) >= big) {
942  big=dum;
943  imax=i;
944  }
945  }
946 
947  //Do we need to interchange rows?
948  if (j != imax) {
949  //Yes, do so...
950  for (k=1;k<=n;k++) {
951  dum=A(imax-1,k-1);
952  A(imax-1,k-1)=A(j-1,k-1); // Tadej: imax-1,k looks wrong
953  A(j-1,k-1)=dum;
954  }
955  //...and change the parity of d.
956  d = -d;
957  vv[imax-1]=vv[j-1]; //Also interchange the scale factor.
958  }
959  indx[j-1]=imax;
960 
961  //If the pivot element is zero the matrix is singular (at least to the precision of the
962  //algorithm). For some applications on singular matrices, it is desirable to substitute
963  //TINY for zero.
964  if (A(j-1,j-1) == 0.0) A(j-1,j-1)=1e-20;
965 
966  //Now, finally, divide by the pivot element.
967  if (j != n) {
968  dum=1.0/(A(j-1,j-1));
969  for (i=j+1;i<=n;i++) A(i-1,j-1) *= dum;
970  }
971  } //Go back for the next column in the reduction.
972 }
#define Assert(Cond)
Definition: bd.h:251
static void nrerror(const TStr &error_text)
Definition: linalg.cpp:655
TSizeTy GetRows() const
Definition: ds.h:2252
static double Abs(const double &Flt)
Definition: dt.h:1427
void Reserve(const TSizeTy &_MxVals)
Reserves enough memory for the vector to store _MxVals elements.
Definition: ds.h:543
TSizeTy GetCols() const
Definition: ds.h:2253
void TNumericalStuff::LUSolve ( const TFltVV A,
const TIntV indx,
TFltV b 
)
static

Definition at line 974 of file linalg.cpp.

974  {
975  Assert(A.GetRows() == A.GetCols());
976  int n = A.GetRows();
977  int i,ii=0,ip,j;
978  double sum;
979  for (i=1;i<=n;i++) {
980  ip=indx[i-1];
981  sum=b[ip-1];
982  b[ip-1]=b[i-1];
983  if (ii)
984  for (j=ii;j<=i-1;j++) sum -= A(i-1,j-1)*b[j-1];
985  else if (sum) ii=i;b[i-1]=sum;
986  }
987  for (i=n;i>=1;i--) {
988  sum=b[i-1];
989  for (j=i+1;j<=n;j++) sum -= A(i-1,j-1)*b[j-1];
990  b[i-1]=sum/A(i-1,i-1);
991  }
992 }
#define Assert(Cond)
Definition: bd.h:251
TSizeTy GetRows() const
Definition: ds.h:2252
TSizeTy GetCols() const
Definition: ds.h:2253
void TNumericalStuff::nrerror ( const TStr error_text)
staticprivate

Definition at line 655 of file linalg.cpp.

655  {
656  printf("NR_ERROR: %s", error_text.CStr());
657  throw new TNSException(error_text);
658 }
char * CStr()
Definition: dt.h:476
double TNumericalStuff::pythag ( double  a,
double  b 
)
staticprivate

Definition at line 660 of file linalg.cpp.

660  {
661  double absa = fabs(a), absb = fabs(b);
662  if (absa > absb)
663  return absa*sqrt(1.0+sqr(absb/absa));
664  else
665  return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+sqr(absa/absb)));
666 }
static double sqr(double a)
Definition: linalg.cpp:647
double TNumericalStuff::sign ( double  a,
double  b 
)
staticprivate

Definition at line 651 of file linalg.cpp.

651  {
652  return b >= 0.0 ? fabs(a) : -fabs(a);
653 }
void TNumericalStuff::SolveLinearSystem ( TFltVV A,
const TFltV b,
TFltV x 
)
static

Definition at line 994 of file linalg.cpp.

994  {
995  TIntV indx; double d;
996  LUDecomposition(A, indx, d);
997  x = b;
998  LUSolve(A, indx, x);
999 }
static void LUSolve(const TFltVV &A, const TIntV &indx, TFltV &b)
Definition: linalg.cpp:974
static void LUDecomposition(TFltVV &A, TIntV &indx, double &d)
Definition: linalg.cpp:909
void TNumericalStuff::SolveSymetricSystem ( TFltVV A,
const TFltV b,
TFltV x 
)
static

Definition at line 837 of file linalg.cpp.

837  {
838  IAssert(A.GetRows() == A.GetCols());
839  TFltV p; CholeskyDecomposition(A, p);
840  CholeskySolve(A, p, b, x);
841 }
#define IAssert(Cond)
Definition: bd.h:262
static void CholeskyDecomposition(TFltVV &A, TFltV &p)
Definition: linalg.cpp:797
TSizeTy GetRows() const
Definition: ds.h:2252
static void CholeskySolve(const TFltVV &A, const TFltV &p, const TFltV &b, TFltV &x)
Definition: linalg.cpp:815
TSizeTy GetCols() const
Definition: ds.h:2253
double TNumericalStuff::sqr ( double  a)
staticprivate

Definition at line 647 of file linalg.cpp.

647  {
648  return a == 0.0 ? 0.0 : a*a;
649 }
void TNumericalStuff::SymetricToTridiag ( TFltVV a,
int  n,
TFltV d,
TFltV e 
)
static

Definition at line 668 of file linalg.cpp.

668  {
669  int l,k,j,i;
670  double scale,hh,h,g,f;
671  for (i=n;i>=2;i--) {
672  l=i-1;
673  h=scale=0.0;
674  if (l > 1) {
675  for (k=1;k<=l;k++)
676  scale += fabs(a(i-1,k-1).Val);
677  if (scale == 0.0) //Skip transformation.
678  e[i]=a(i-1,l-1);
679  else {
680  for (k=1;k<=l;k++) {
681  a(i-1,k-1) /= scale; //Use scaled a's for transformation.
682  h += a(i-1,k-1)*a(i-1,k-1);
683  }
684  f=a(i-1,l-1);
685  g=(f >= 0.0 ? -sqrt(h) : sqrt(h));
686  IAssertR(_isnan(g) == 0, TFlt::GetStr(h));
687  e[i]=scale*g;
688  h -= f*g; //Now h is equation (11.2.4).
689  a(i-1,l-1)=f-g; //Store u in the ith row of a.
690  f=0.0;
691  for (j=1;j<=l;j++) {
692  // Next statement can be omitted if eigenvectors not wanted
693  a(j-1,i-1)=a(i-1,j-1)/h; //Store u=H in ith column of a.
694  g=0.0; //Form an element of A  u in g.
695  for (k=1;k<=j;k++)
696  g += a(j-1,k-1)*a(i-1,k-1);
697  for (k=j+1;k<=l;k++)
698  g += a(k-1,j-1)*a(i-1,k-1);
699  e[j]=g/h; //Form element of p in temporarily unused element of e.
700  f += e[j]*a(i-1,j-1);
701  }
702  hh=f/(h+h); //Form K, equation (11.2.11).
703  for (j=1;j<=l;j++) { //Form q and store in e overwriting p.
704  f=a(i-1,j-1);
705  e[j]=g=e[j]-hh*f;
706  for (k=1;k<=j;k++) { //Reduce a, equation (11.2.13).
707  a(j-1,k-1) -= (f*e[k]+g*a(i-1,k-1));
708  Assert(_isnan(static_cast<double>(a(j-1,k-1))) == 0);
709  }
710  }
711  }
712  } else
713  e[i]=a(i-1,l-1);
714  d[i]=h;
715  }
716  // Next statement can be omitted if eigenvectors not wanted
717  d[1]=0.0;
718  e[1]=0.0;
719  // Contents of this loop can be omitted if eigenvectors not
720  // wanted except for statement d[i]=a[i][i];
721  for (i=1;i<=n;i++) { //Begin accumulation of transformationmatrices.
722  l=i-1;
723  if (d[i]) { //This block skipped when i=1.
724  for (j=1;j<=l;j++) {
725  g=0.0;
726  for (k=1;k<=l;k++) //Use u and u=H stored in a to form PQ.
727  g += a(i-1,k-1)*a(k-1,j-1);
728  for (k=1;k<=l;k++) {
729  a(k-1,j-1) -= g*a(k-1,i-1);
730  Assert(_isnan(static_cast<double>(a(k-1,j-1))) == 0);
731  }
732  }
733  }
734  d[i]=a(i-1,i-1); //This statement remains.
735  a(i-1,i-1)=1.0; //Reset row and column of a to identity matrix for next iteration.
736  for (j=1;j<=l;j++) a(j-1,i-1)=a(i-1,j-1)=0.0;
737  }
738 }
#define IAssertR(Cond, Reason)
Definition: bd.h:265
#define Assert(Cond)
Definition: bd.h:251
TStr GetStr() const
Definition: dt.h:1459

The documentation for this class was generated from the following files: