Page MenuHomeHEPForge
Files F8308858

ScannerSAux.cpp
No OneTemporary
Actions

ScannerSAux.cpp
View Options

#include "ScannerSAux.h"
//#include "ScannerSModelclasses.h"
#ifdef VERBOSE
#include "ScannerSVerbose.h"
#endif
using namespace std;
int GaussUpTriang(std::vector< std::vector<double> > & M){
////Initialise counter of horizontal continuations in gaussian elimination
int hcont=0;
///Initialise counter to scan lines vertically
int j=0;
for(int a=0; (a<M[0].size() &&j<M.size()); ++a){
//logical switch to keep finding a line with non-zero first value
bool Continue = true;
//Find a line starting with a nonzero value in this column
do{
if(M[j][a]!=0){
//Line found. Divide all elements by this number so that the leftmost element becomes 1.
for(int a1=M[0].size()-1;a1>a-1; a1--){
M[j][a1]/=M[j][a];
}
//Push this line up
M[j].swap(M[a-hcont]);
//swap(M[j],M[a-hcont]);
//Stop search
Continue=false;
}
//Otherwise continue search
j++;
}while(Continue && j<M.size());
//set j to the vertical position where the line was stored
j=a-hcont;
// If a line was found then perform gaussian elimination on the lines below it
if(Continue==false){
//Loop over lines starting at the one just below
for(int j1=j+1;j1<M.size();j1++){
//Loop over element of the current line to set the first element to zero and recompute the others
for(int a1=M[0].size()-1;a1>a-1;a1--){
//subtract the leading line
double temp=M[j1][a1]-M[j1][a]*M[j][a1];
if(abs(temp)>abs(M[j1][a1])*1e-10)
M[j1][a1]=temp;
else
M[j1][a1]=0e0;
}
}
}
else{
//Otherwise if all elements in that sub-column are zero increase count of horizontal continuations
hcont++;
}
//Compute j and add 1 to follow increase in a for the next loop
j=a-hcont+1;
}
return hcont;
}
void GaussSolveUpTriang(std::vector< std::vector<double> > & M,unsigned int & NIC,std::vector<unsigned int> & order, int hcont){
// Find position of variables to be eliminated from the system
unsigned int count=0;
//Scan over (partially) upper triangular section
for(int j=0; (j< M[0].size()-hcont && j<M.size());j++){
//Scan line from left to right until we hit a 1.
for(int a=j; a< M[0].size(); a++){
if(M[j][a]==1){
//Store position of this 1.
order[count]=a;
a=M[0].size();
count++;
}
}
}
NIC=count;
//Store positions of remaining couplings at the end of order vector
count=0;
int acount=0;
for(int a=0;a<M[0].size(); a++){
if(a!=order[count]){
order[NIC+acount]=a;
acount++;
}
else{
count++;
}
}
////Loop over upper tringular from bottom to top line by line
int jmax=M[0].size()-hcont-1;
if(jmax>M.size()-1)
jmax=M.size()-1;
for(int j=jmax; j>0;j--){
////Loop over current line from right to left to create a zero at the corresponding position
for(int j1=jmax; j1>j-1; j1--){
////eliminate element using the correct line below and calculate line from the correct position a to the right
for(int a=M[0].size()-1; a> order[j1]-1; a--){
M[j-1][a]-=M[j-1][order[j1]]*M[j1][a];
}
}
}
}
void RedBydV(vector<vector<double> > & dV,deriv2V & d2V, unsigned int NIC, std::vector<unsigned int> & orderdV){
for(size_t a2=NIC; a2!=d2V[0][0].size(); a2++){
d2V.Mag[a2]=0e0;
//Reduce second order conditions by eliminating the first order coefficients
for(size_t i=0; i!=d2V.size(); i++){
for(size_t j=i; j!=d2V.size(); j++){
double sum=0e0;
for(size_t a1=0; a1!=NIC; a1++){
sum+=d2V[i][j][orderdV[a1]]*dV[a1][orderdV[a2]];
}
double temp=d2V[i][j][orderdV[a2]]-sum;
if(abs(temp)>abs(d2V[i][j][orderdV[a2]])*1e-10)
d2V[i][j][orderdV[a2]]-=sum;
else
d2V[i][j][orderdV[a2]]=0e0;
d2V[j][i][orderdV[a2]]=d2V[i][j][orderdV[a2]];
// Store value so that it serves as a scale later...
if(abs(d2V[i][j][orderdV[a2]])>d2V.Mag[a2])
d2V.Mag[a2]=abs(d2V[i][j][orderdV[a2]]);
}
}
}
}
void FindCandidates(deriv2V & d2V,double * & TempM, double * & TempOut,std::vector<unsigned int> & orderdV,unsigned int & a2ref,unsigned int & kflat,unsigned int NIC){
bool kbool=true;
for(size_t a2=NIC; a2!=d2V[0][0].size() && kbool;++a2){
// Analysing flat direction candidates
for(size_t i=0; i!=d2V.size(); ++i)
for(size_t j=0; j!= d2V.size(); ++j)
TempM[d2V.size()*i+j]=d2V[i][j][orderdV[a2]];
/// Finding eigenvectors and eigenvalues of reference matrix
EigenSystemRsym(d2V.size(),TempM,TempOut);
for(unsigned int j=0;j!=d2V.size();++j){
//found a non-zero eigenvalue
if(abs(TempOut[j])>1e-10*abs(TempOut[d2V.size()-1])&& kbool){
kflat=j;
kbool=false;
}
}
a2ref=a2;
}
}
void TestCandidates(deriv2V & d2V, double * & TempOut,std::vector<unsigned int> & orderdV,std::vector<unsigned int> & Gpos,unsigned int & kminus, unsigned int & kminus2, unsigned int kflat,unsigned int a2ref,unsigned int NIC){
//Check other couplings by computing matrix elements in basis of the first coupling
vector<double> Vtest(d2V.size()); // Auxiliary vector to test matrix elements
for(size_t a2=a2ref+1; a2!=d2V[0][0].size(); ++a2){
for(size_t k=0;k!=kflat; ++k){
//Check if already excluded (does nothing in the first call of the loop since kminus=0 or if kminus remains zero)
bool skip=false;
for(size_t k1=0; k1!=kminus; ++k1){
if(k==Gpos[k1])
skip=true;
}
//Compute action of d2V[i][j][] on k candidate
for(size_t i=0;i!=d2V.size() && (not skip);++i){
Vtest[i]=0e0;
for(size_t j=0; j!=d2V.size(); ++j)
Vtest[i]+=d2V[i][j][orderdV[a2]]*TempOut[d2V.size()*(k+1)+j];
//Test if still an eigenvector with zero eigenvalue
if(abs(Vtest[i])>1e-10*d2V.Mag[a2]){
//If component not consistent, then store position of vector that failed...
Gpos[kminus]=k;
// and increase count of candidates which fail...
kminus++;
i=d2V.size()-1; // skip checks of remaining components since the vector has failed
}
}
}
}
// Directions which are not flat and may be curved eigen-drections
unsigned int kleft=kminus+d2V.size()-kflat;
kminus2=kleft;//kminus2 will contain the number of directions that will mix. So in the worst case everything left will mix
if(kleft>1){
double * Mtest;
double * MtestOut;
vector<double> tempvec(d2V.size(),0);
vector< vector<double> > Basis(kleft,tempvec),BasisStore(kleft,tempvec); //Basis composed of eigenvectors of the reference matrix, and Store Space for basis vectors as they are accepted or rejected
Mtest= new double [kleft*kleft]; //Space for reference matrix being checked
MtestOut= new double [(kleft+1)*kleft]; //Eigenvalues and eigenvectors of the reference matrix
//Initialise basis using the first reference matrix that was used when finding candidates
////First place the rejected flat directions of the reference matrix
for(size_t k=0; k!=kminus; ++k)
for(size_t j=0; j!=d2V.size(); ++j)
Basis[k][j] = TempOut[d2V.size()*(Gpos[k]+1)+j];
////Then place the curved directions
for(size_t k=kflat; k!=d2V.size(); ++k)
for(size_t j=0; j!=d2V.size(); ++j)
Basis[kminus+k-kflat][j] = TempOut[d2V.size()*(k+1)+j];
unsigned int a2ref2=a2ref; //Integer to loop over other matrices to check their eigenvectors
unsigned int kaccept=0; //Counter of the curved directions that are accepted
unsigned int dimN=kleft; //Dimension of the subspace under inspection
while(a2ref2<d2V[0][0].size() && dimN>0){//Use matrix a2ref2 as reference if there is a subspace to check
//Initialise the reference matrix a2ref2 in the basis of the subspace being checked
for(size_t k1=0; k1!= dimN; ++k1){
for(size_t k2=0; k2!=dimN; ++k2){
Mtest[k1*dimN+k2] = 0e0;
for(size_t i=0; i!=d2V.size(); ++i){
double TempSum=0e0;
for(size_t j=0; j!=d2V.size(); ++j)
TempSum+=d2V[i][j][orderdV[a2ref2]]*Basis[k2][j];
Mtest[k1*dimN+k2]+= Basis[k1][i]*TempSum;
}
}
}
//Compute the eigenvalues of the reference matrix in this basis
EigenSystemRsym(dimN,Mtest,MtestOut);
size_t kfail=0; //Integer to count vectors that are rejected
//Loop over all vectors to check if any can be a curved eigendirection
for(size_t k=0;k!=dimN;++k){
// Expand back the eigenvector in the original basis
for(size_t i=0; i!=d2V.size(); ++i){
tempvec[i]=0e0;
for(size_t k2=0;k2!=dimN;++k2)
tempvec[i]+=MtestOut[(k+1)*dimN+k2]*Basis[k2][i];
if(abs(tempvec[i])<1e-10)
tempvec[i]=0e0;//Chop numerical errors
}
//Check this eigenvector if the eigenvalue is non-zero, to see if it can be a curved eigendirection
if(MtestOut[k]!=0){
///test other matrices against such vector
bool testcurv=false; //logical switch to check if the direction has been accepted
for(size_t a2=a2ref;a2!=d2V[0][0].size();++a2){
//loop over all other matrices
if(a2!=a2ref2){
//test vector against matrix
testcurv=TestCurved(a2,NIC,d2V,orderdV,tempvec);
if(not testcurv){
//Direction failed to be accepted....
for(size_t i=0;i!=d2V.size();++i){
if(abs(tempvec[i])>1e-10)
BasisStore[kfail][i]=tempvec[i]; //Store basis in the first range from k=0 with increasing kfail as failed directions are spotted
else
BasisStore[kfail][i]=0e0; //Chop numerical errors (note the condition above is good because the eigenvectors are normalised)
}
kfail++; //Add to the counter of vectors that have been excluded
a2=d2V[0][0].size()-1; // avoid check of other matrices since the direction is already excluded
}
}
}
// Case when the direction was accepted in the previous loop
if(testcurv){
kaccept++; //increase counter of accepted curved directions
for(size_t i=0;i!=d2V.size();++i){//And store vector in the top end of the k range (rejected directions are stored down up, and accepted directions, top down)
if(abs(tempvec[i])>1e-10)
BasisStore[kleft-kaccept][i]=tempvec[i];
else
BasisStore[kleft-kaccept][i]=0e0; //Chop numerical errors
}
}
}
else{//Otherwise its a "null" eigendirection so we cannot accept it as a curved direction or reject it (because it may be an eigendirection of one or more of the other matrices with non-zero eigenvalue)
for(size_t i=0;i!=d2V.size();++i){
//Store direction that may or may not be curved
if(abs(tempvec[i])>1e-10)
BasisStore[kfail][i]=tempvec[i];
else
BasisStore[kfail][i]=0e0;//Chop numerical errors
}
kfail++; //Icrease count od directions of the new subspace to be checked in the next round of a2ref2
}
}
dimN=kleft-kaccept;// Dimension of new subspace to be checked
//Set new basis for (possibly) smaller subspace using the stored basis vectors
for(size_t k=0; k!=dimN; ++k)
for(size_t j=0; j!=d2V.size(); ++j)
Basis[k][j] = BasisStore[k][j];
a2ref2++;
}
//The new basis has been defined with the directions that will mix in the lower sector of BasisStore and the curved eigendirections in the upper sector of the array (in k)
//// Replace the vectors in TempOut which were candidate flat directions that failed the tests by the first few vectors in BasisStore
for(size_t k=0;k!=kminus;++k)
for(size_t j=0; j!=d2V.size(); ++j)
TempOut[d2V.size()*(Gpos[k]+1)+j]=BasisStore[k][j];
//// Replace the remaining vectors which were candidate curved directions by the remaining vectors in BasisStore
for(size_t k=kflat;k!=d2V.size();++k)
for(size_t j=0; j!=d2V.size(); ++j)
TempOut[d2V.size()*(k+1)+j]=BasisStore[kminus+k-kflat][j];
delete [] Mtest;
delete [] MtestOut;
//Set kminus2 to contain the number of directions that mix
kminus2=kleft-kaccept;
}
else{
kminus2=0;
}
}
bool TestCurved(size_t a2,size_t NIC,deriv2V & d2V, std::vector<unsigned int> & orderdV,std::vector<double> & vec){
vector<double> Vtest(d2V.size(),0); // Auxiliary vector to test matrix elements
bool result=true;
//Compute estimate for eigenvalue
double lambdaEst=0e0;
for(size_t i=0;i!=d2V.size();++i){
Vtest[i]=0e0;
for(size_t j=0; j!=d2V.size(); ++j)
Vtest[i]+=d2V[i][j][orderdV[a2]]*vec[j];
lambdaEst+=Vtest[i]*vec[i]; //Assuming vec[i] is normalised, if we are dealing with an eigendirection then the dot product give the eigenvalue
}
//Compute action of d2V[i][j][] on k candidate
for(size_t i=0;i!=d2V.size();++i){
Vtest[i]=0e0;
for(int j=0; j!=d2V.size();++j)
Vtest[i]+=d2V[i][j][orderdV[a2]]*vec[j];
///Use some scales to assess if candidate vector component is non-zero (otherwise we cannot rule out it is an eigenvector from this component)
if(abs(lambdaEst)>1e-10*d2V.Mag[a2] && abs(vec[i])>1e-10*Vtest[i]/lambdaEst){
double lambda=Vtest[i]/vec[i]; //lambda obtained from this component
//Check if the current component is consistent with the eigenvalue estimate
if(abs(lambdaEst-lambda)>1e-10*d2V.Mag[a2]){
result=false;//If no, i.e. estimate differs too much from the current component then reject
i=d2V.size()-1;// avoid continuing the computation of other matrix elements
}
}
else if(abs(Vtest[i])>1e-10*d2V.Mag[a2]){//In this case lambda should be zero so it is inconsistent to have a non-zero value for any of the components of Vtest. Then reject candidate.
result=false;
i=d2V.size()-1; // cancel further checks of matrix elements
}
}
return result;
}
void ConvertEigen(deriv2V & d2V,deriv2V & d2Vnew,double * & TempM, double * & TempOut,std::vector<unsigned int> & orderdV,std::vector<unsigned int> & Gpos,unsigned int kminus, unsigned int kminus2,unsigned int kflat,unsigned int a2ref,unsigned int NIC){
///Converting system to basis with a priori diagonal directions. The basis will be composed first by the vectors which are not eigenvectors and then the a priori eigen-vectors
//// If there are zero matrices before the reference matrix, set them to zero
for(size_t a2=NIC; a2!=a2ref;++a2)
for(size_t i=0; i!=d2V.size(); ++i)
for(size_t j=0; j!=d2V.size(); j++)
d2Vnew[i][j][orderdV[a2]]=0e0;
//// Counter of non-eigen directions
unsigned int kNo=0,kYesF=0,kYesC=0;
/// Place non-eigen directions first in the matrix
for(size_t i=0; i!=d2V.size(); ++i){
//Test if the ith direction is an eigen direction
bool accept=true;
bool flat=false;
if(i<kflat)
flat=true;
//Loop over non-eigen directions
if(kminus>kminus2){//Case when the number of rejected flat directions is larger than the number of directions that will mix which means that all mixing directions are part of this set and other directions are curved eigendirections
for(size_t k1=0;k1!=kminus;++k1){
//reject ith direction in case it's on the list Gpos
if(Gpos[k1]==i)
accept=false;
else
flat=true;
}
}
else{//Complementary case when some of the mixing directions are in the set counted by kminus and others are from the initial candidate curved eigendirections
for(size_t k1=0;k1!=kminus;++k1){
//reject ith direction in case it's on the list Gpos
if(Gpos[k1]==i)
accept=false;
else
flat=true;
}
for(int k1=kflat;k1!=kflat+kminus2-kminus;++k1){
//reject ith direction in case it has been rejected in the list of potential curved eigen-directions
if(k1==i)
accept=false;
}
}
// If rejected then place the corresponding eigenvalue in the kth position and, increase the count by 1 and place the eigenvector in TempM for simplicity
if(not accept){
for(size_t s=0;s!=d2V.size();++s){
TempM[d2V.size()*kNo+s]=TempOut[d2V.size()*(i+1)+s];
}
kNo++;
}
else{//Otherwise place the eigenvalues in reverse order in the remaining diagonals and put eigenvectors in TempM in the same reversed order
if(flat){
for(size_t s=0;s!=d2V.size();++s)
TempM[d2V.size()*(d2V.size()-1-kYesF)+s]=TempOut[d2V.size()*(i+1)+s];
kYesF++;
}
else{
for(size_t s=0;s!=d2V.size();++s)
TempM[d2V.size()*(d2V.size()-1-(kflat-kminus)-kYesC)+s]=TempOut[d2V.size()*(i+1)+s];
kYesC++;
}
}
}
///////////////////////////////////////////////////////////////
for(size_t a2=a2ref;a2!=d2V[0][0].size();++a2){
//// Secondly compute the matrix elements all matrices
vector<double> Vtest(d2V.size(),0); // Auxiliary vector to test matrix elements
for(size_t i=0; i!=d2V.size(); ++i){
//Compute action on eigenvector
for(size_t k=0;k!=d2V.size();++k){
Vtest[k]=0e0;
for(size_t j=0; j!=d2V.size(); ++j){
Vtest[k]+=d2V[k][j][orderdV[a2]]*TempM[d2V.size()*i+j];
}
}
//Compute matrix elements
for(size_t k=i;k!=d2V.size();++k){
double sum=0e0;
for(size_t l=0; l!=d2V.size();++l){
sum+=TempM[d2V.size()*k+l]*Vtest[l];
}
d2Vnew[k][i][orderdV[a2]]=sum;
d2Vnew[i][k][orderdV[a2]]=sum;
}
}
}
//Chop rounding errors for eigen-directions
for(size_t a2=NIC;a2!=d2V[0][0].size();++a2){
for(size_t i= kminus+kminus2; i!=d2V.size(); ++i){
for(size_t k=0;k!=d2V.size();++k){
if(k!=i){
d2Vnew[k][i][orderdV[a2]]=0e0;
d2Vnew[i][k][orderdV[a2]]=0e0;
}
else if(abs(d2Vnew[k][i][orderdV[a2]])<1e-10*d2V.Mag[a2]){
d2Vnew[k][i][orderdV[a2]]=0e0;
}
}
}
}
}
void GenerateMrot(deriv2V & d2Vnew,std::vector< std::vector<double> > & Mrot,std::vector<unsigned int> & orderdV,unsigned int kminus2,unsigned int NIC,RandGen & r){
//Generating random rotation matrix
double *Mdata, *Mdata2, *Mdata3,*bdata;
Mdata = new double [(kminus2)*(kminus2)];
Mdata2 = new double [(kminus2)*(kminus2)];
Mdata3 = new double [(kminus2)*(kminus2)];
bdata = new double [kminus2];
for(size_t i=0;i!=kminus2; ++i){
for(size_t j=0;j!=kminus2; ++j){
//Initialise to a normal distributed random matrix
Mrot[i][j]=2*r(0,1)-1e0;
while(Mrot[i][j]==0){
// Avoid the rare case that we hit a zero exactly (this would cause numerical problems in the next step)
Mrot[i][j]=2*r(0,1)-1e0;
}
if(Mrot[i][j]<0)
Mrot[i][j]=-sqrt(2e0*log(1e0/abs(Mrot[i][j])));
else
Mrot[i][j]=sqrt(2e0*log(1e0/abs(Mrot[i][j])));
Mdata[(kminus2)*i+j]=Mrot[i][j];
}
}
#ifdef VERBOSE
clog << "Printing out Gaussian matrix used to obtain Mmixing..." << endl;
PrintM1(Mrot);
#endif
size_t tempdim=kminus2;
if(tempdim==0)
tempdim=1;
gsl_matrix_view m = gsl_matrix_view_array(Mdata,tempdim,tempdim);
// exit(0);
gsl_vector_view b = gsl_vector_view_array(bdata,tempdim);
gsl_linalg_QR_decomp (&m.matrix, &b.vector);
gsl_matrix_view Q = gsl_matrix_view_array(Mdata2,tempdim,tempdim);
gsl_matrix_view R = gsl_matrix_view_array(Mdata3,tempdim,tempdim);
////Define a temporary matrix d2Vtemp
deriv2V d2Vtemp=d2Vnew;
if(kminus2!=0){
gsl_linalg_QR_unpack(&m.matrix,&b.vector, &Q.matrix, &R.matrix);
for(size_t i=0;i!=kminus2; ++i)
for(int j=0;j!=kminus2; ++j)
Mrot[i][j]=gsl_matrix_get(&Q.matrix,i,j);
if(2*(kminus2/2)==kminus2)
Mrot[0].swap(Mrot[1]);
#ifdef VERBOSE
clog<< "Printing out the sub-mixing matrix..."<< endl;
PrintM1(Mrot);
#endif
//Rotate d2VnewMatrix
for(size_t a2=NIC; a2!=d2Vnew[0][0].size(); ++a2){
for(size_t i=0;i!=kminus2;++i){
for(size_t j=0;j!=kminus2;++j){
double sum=0e0;
for(size_t k=0;k!=kminus2;++k)
sum+=d2Vnew[i][k][orderdV[a2]]*Mrot[j][k];
d2Vtemp[i][j][orderdV[a2]]=sum;
}
}
}
for(size_t a2=NIC; a2!=d2Vnew[0][0].size(); ++a2){
for(size_t i=0;i!=kminus2;++i){
for(size_t j=0;j!=kminus2;++j){
double sum=0e0;
for(size_t k=0;k!=kminus2;++k){
sum+=Mrot[i][k]*d2Vtemp[k][j][orderdV[a2]];
}
d2Vnew[i][j][orderdV[a2]]=sum;
}
}
}
}
#ifdef VERBOSE
clog << endl<< "****** d2Vnew independent matrices in rotated basis ******" << endl;
PrintM22(d2Vnew,NIC,orderdV);
#endif
delete [] Mdata;
delete [] Mdata2;
delete [] Mdata3;
delete [] bdata;
}
//////////////////
void PrintNflatCandidates(deriv2V & d2V,double * & TempM,double * & TempOut,std::vector<unsigned int> & orderdV,unsigned int NIC){
for(size_t a2=NIC; a2!= d2V[0][0].size(); ++a2){
for(size_t i=0; i!=d2V.size(); ++i)
for(size_t j=0; j!=d2V.size();++j)
TempM[d2V.size()*i+j]=d2V[i][j][orderdV[a2]];
EigenSystemRsym(d2V.size(),TempM,TempOut);
clog << endl << "----Eigenvalues for matrix a ="<< orderdV[a2]<< endl;
for(size_t i=0; i!=d2V.size(); ++i){
clog << TempOut[i]<<"\t";
}
clog << endl << "----Eigenvectors for matrix a ="<< orderdV[a2]<< endl;
for(size_t i=0; i!=d2V.size(); ++i){
for(size_t j=0;j!=d2V.size(); ++j)
clog << TempOut[(i+1)*d2V.size()+j]<<"\t";
clog << endl;
}
if(abs(TempOut[d2V.size()-1])>1e-10){
size_t k;
for(size_t j=0;j!=d2V.size();++j){
if(abs(TempOut[j])>1e-10*abs(TempOut[d2V.size()-1])){
k=j;
j=d2V.size()-1;
}
}
clog << endl << k << " Candidate flat directions "<< endl;
}
else{
clog << d2V.size() << " Candidate flat directions "<< endl;
}
}
clog << endl;
}
///////////////////////////////////////////////////////////////////////
///////////////// GENERIC MATHEMATICAL ALGORITHMS /////////////////////
///////////////////////////////////////////////////////////////////////
///Solve the eigenvalue problem (adapted from gsl examples)
void EigenSystemRsym(int N,double * & DataIn, double * & DataOut){
int i,j;
gsl_matrix_view m1 = gsl_matrix_view_array (DataIn, N, N);
gsl_vector *eval = gsl_vector_alloc(N);
gsl_matrix *evec = gsl_matrix_alloc(N, N);
gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc(N);
gsl_eigen_symmv (&m1.matrix, eval, evec, w);
gsl_eigen_symmv_free (w);
gsl_eigen_symmv_sort (eval, evec,GSL_EIGEN_SORT_ABS_ASC);
for (i = 0; i < N; i++)
{
double eval_i = gsl_vector_get (eval, i);
DataOut[i]= eval_i;
gsl_vector_view evec_i = gsl_matrix_column (evec, i);
for (j = 0; j < N; j++)
{
double z = gsl_vector_get(&evec_i.vector, j);
DataOut[N+N*i+j]=z;
}
}
gsl_vector_free(eval);
gsl_matrix_free(evec);
}
///Solve the eigenvalue problem but no sorting and no vectors on output!!! (adapted from gsl examples)
void EigenSystemRsymNoSort(int N,double * & DataIn, double * & DataOut){
int i,j;
gsl_matrix_view m1 = gsl_matrix_view_array (DataIn, N, N);
gsl_vector *eval = gsl_vector_alloc(N);
gsl_matrix *evec = gsl_matrix_alloc(N, N);
gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc(N);
gsl_eigen_symmv (&m1.matrix, eval, evec, w);
gsl_eigen_symmv_free (w);
//gsl_eigen_symmv_sort (eval, evec,GSL_EIGEN_SORT_ABS_ASC);
for (i = 0; i < N; i++)
{
double eval_i = gsl_vector_get (eval, i);
DataOut[i]= eval_i;
}
gsl_vector_free(eval);
gsl_matrix_free(evec);
}
double factorial(double n){
double temp=1;
for(size_t i=1;i<n;++i){
temp*=i+1;
}
return temp;
}
double Integrator(double A, double (*Integrand)(double, void *))
{
gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);
double result, error;
//A_params params = {A};
// Assign gsl function to integrate
gsl_function F;
F.function = Integrand;
F.params = &A;
// perform the integration using an adaptive method which can deal with singularities.
gsl_integration_qags(&F, 0, 1, 0, 0.000001, 1000,w, &result, &error);
// Free memory
gsl_integration_workspace_free (w);
return result;
}
double Interp(double x, double * table){
unsigned int dim=(unsigned int)table[0];
double xmin=table[1];
double xmax=table[2];
double dix=(x-xmin)*(dim-1)/(xmax-xmin);
unsigned int ix=(unsigned int)dix;
if(dix< 0 || dix>dim-1)
cerr<< "ScannerSAux.cpp: double Interp(double x, double * table). Error, interpolation point out of bounds. Check your code!!!"<< endl;
return table[ix+3]+(table[ix+4]-table[ix+3])*(dix-ix);
}

File Metadata

Mime Type
text/x-c
Expires
Sat, Dec 21, 1:24 PM (1 d, 1 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4008927
Default Alt Text
ScannerSAux.cpp (23 KB)

Event Timeline