Page MenuHomeHEPForge
Files F8309803

ScannerSAux.cpp
No OneTemporary
Actions

ScannerSAux.cpp
View Options

#include "ScannerSAux.h"
//#include "ScannerSModelclasses.h"
#ifdef VERBOSE
#include "ScannerSVerbose.h"
#endif
const double epsilon=1e-10;
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])*epsilon)
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]])*epsilon)
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 FindAprioriEigenDirs(deriv2V & d2V, std::vector< std::vector<double> > & Mbasis,double * Agen,double * Bgen,double * MtestOut,std::vector<unsigned int> & orderdV, std::vector<int> & degen,vector<int> & BlockInfo,vector<unsigned int> & BlockDims, unsigned int & kflat,unsigned int & kcurved,unsigned int NIC){
RandGen r;
degen.resize(0); //empty the vector which will hold the multiplicities of the degenerate curved directions. This will be in an inverted order due to the storage order of the eigendirections as they are found.
BlockDims.resize(0);
//Creating auxiliary generic matrix...
unsigned int adim=d2V[0][0].size()-NIC; //Number of independent couplings
vector<double> randvec(adim,0); // generic random coupling vector
for(size_t a=0;a<adim;++a){
randvec[a]=r(-1,1); //creating random couplings...
if(d2V.Mag[a+NIC]!=0)
randvec[a]/=d2V.Mag[a+NIC]; //Using typical scale of matrix elements for this matrix to re-scale random number
}
////filling generic matrix through the linear combination...
for(size_t k=0; k!=d2V.size(); ++k){
Agen[k*d2V.size()+k] = 0e0; //Set diagonal element to zero
for(size_t a=0;a!=adim;++a)//filling diagonal element through linear combination
Agen[k*d2V.size()+k] +=d2V[k][k][orderdV[a+NIC]]*randvec[a];
for(size_t j=k+1; j!=d2V.size(); ++j){
Agen[j*d2V.size()+k] = 0e0; //Set off-diagonal element to zero
for(size_t a=0;a!=adim;++a)
Agen[j*d2V.size()+k]+=d2V[j][k][orderdV[a+NIC]]*randvec[a]; //filling off-diagonal element
Agen[k*d2V.size()+j] = Agen[j*d2V.size()+k] ; //Use that matrices are symmetric to get transposed term
}
}
//Compute the eigenvalues of the reference matrix in this basis
EigenSystemRsym(d2V.size(),Agen,MtestOut);
double ScaleEigen=1;//Variable to hold a scale for the eigenvalues, to test for null eigenvalues.
if(abs(MtestOut[d2V.size()-1])!=0)
ScaleEigen=MtestOut[d2V.size()-1];//If not all eigenvalues are zero, use largest one to set the scale.
for(size_t dummy=0;dummy<d2V.size();++dummy){//Loop to detect all flat directions
if(abs(MtestOut[kflat]/ScaleEigen)<epsilon){//Flat direction detected!
kflat++; //increase counter of accepted flat directions
for(size_t j=0;j!=d2V.size();++j){//And store vector at the end of the k range (accepted flat directions are stored from the end of Mbasis as they are found)
if(abs(MtestOut[d2V.size()*kflat+j])>epsilon)//Non-zero components
Mbasis[d2V.size()-kflat][j]=MtestOut[d2V.size()*kflat+j];
else
Mbasis[d2V.size()-kflat][j]=0e0; //Chop numerical errors
}
}
else
break;//Otherwise, since eigenvalues are stored with increasing magnitude, stop the search loop
}
for(size_t a=0;a<adim;++a){
randvec[a]=r(-1,1); //creating random couplings again for second auxiliary matrix Bgen...
if(d2V.Mag[a+NIC]!=0)
randvec[a]/=d2V.Mag[a+NIC]; //Using typical scale of matrix elements for this matrix to re-scale random number
}
// Compute generic Bgen in the basis of the eigenvectors Agen and fill vector with block information
for(size_t i=0;i!=d2V.size();++i)
BlockInfo[i]=0; //Initializing Block info vector to zero to be used to get block directions info
unsigned int labelblock=1; //loop variable to be increased as we move to a new block. 1 means we are in the first "block"
////filling generic matrix through the linear combination (sum over couplings times basis matrices)...
for(size_t k=0; k!=d2V.size()-kflat; ++k){//Loop over all directions
Bgen[k*d2V.size()+k] = 0e0; //Set diagonal element to zero
for(size_t a=0;a!=adim;++a){//filling diagonal element through linear combination
//Compute matrix element in new basis
double Mkk=0;
for(size_t alpha=0;alpha<d2V.size();++alpha){
double Mkkalpha=0;//Variable to hold sum over couplings for component kk
for(size_t beta=0;beta<d2V.size(); ++beta)
Mkkalpha+=d2V[alpha][beta][orderdV[a+NIC]]*MtestOut[d2V.size()*(k+1+kflat)+beta];
Mkk+=MtestOut[d2V.size()*(k+1+kflat)+alpha]*Mkkalpha;
}
Bgen[k*d2V.size()+k] +=Mkk*randvec[a];//Store matrix element
}
if(abs(Bgen[k*d2V.size()+k]/ScaleEigen)<epsilon)
Bgen[k*d2V.size()+k]=0;
bool blockfound=false;//Variable to test if a block was found
for(size_t j=k+1; j!=d2V.size()-kflat; ++j){
Bgen[j*d2V.size()+k] = 0e0; //Set off-diagonal element to zero
for(size_t a=0;a!=adim;++a){
//Compute matrix element in new basis
double Mjk=0;
for(size_t alpha=0;alpha<d2V.size();++alpha){
double Mjkalpha=0;
for(size_t beta=0;beta<d2V.size(); ++beta)
Mjkalpha+=d2V[alpha][beta][orderdV[a+NIC]]*MtestOut[d2V.size()*(k+1+kflat)+beta];
Mjk+=MtestOut[d2V.size()*(j+1+kflat)+alpha]*Mjkalpha;
}
Bgen[j*d2V.size()+k] +=Mjk*randvec[a];//Store matrix element
}
if(abs(Bgen[j*d2V.size()+k]/ScaleEigen)<epsilon)
Bgen[j*d2V.size()+k]=0;//Chop numerical errors
Bgen[k*d2V.size()+j] = Bgen[j*d2V.size()+k] ; //Use that matrices are symmetric to get transposed term
//Finally test if the direction os part of a block (first condition) and if it wasn't added to a block already
if(Bgen[k*d2V.size()+j]!=0 && BlockInfo[k]==0){
blockfound=true;
BlockInfo[j]=labelblock; //Label direction with the number of the block
}
}
if(blockfound && BlockInfo[k]==0){//otherwise eigen-direction or already found
BlockInfo[k]=labelblock;
labelblock++;
}
else if(abs(Bgen[k*d2V.size()+k]/ScaleEigen)>epsilon && not blockfound && BlockInfo[k]==0){//otherwise curved eigen direction or already found
for(int k1=k-1; k1>=0;--k1){
if(abs(Bgen[k*d2V.size()+k]-Bgen[k1*d2V.size()+k1])<epsilon){
BlockInfo[k]=BlockInfo[k1]; //Part of another block so give it the same label
break; //avoid further (un-necessary) checks
}
}
if(BlockInfo[k]==0){
BlockInfo[k]=-labelblock; //curved direction found
labelblock++;
}
}
}
//Organising a priori degenerate curved directions...
for(size_t k=0; k<d2V.size()-kflat;++k){
int tempdegen=1;//variable to count number of degenerate states
if(BlockInfo[k]<0){ //if eigen-directions...
for(size_t j=k+1; j<d2V.size()-kflat;++j){ //Loop over all other directions to find which are degenerate with this one...
if(BlockInfo[j]==BlockInfo[k]){
++kcurved; //increase counter of curved directions
++tempdegen; //increse degeneracy label for this set of curved directions
for(size_t m=0;m!=d2V.size();++m){//Store next degenerate direction...
if(abs(MtestOut[d2V.size()*(j+kflat+1)+m])>epsilon)//Non-zero element
Mbasis[d2V.size()-(kflat+kcurved)][m]=MtestOut[d2V.size()*(j+kflat+1)+m];//place direction from bottom of matrix up on "top" of the flat directions
else
Mbasis[d2V.size()-(kflat+kcurved)][m]=0e0; //Chop numerical errors
}
//Destroy BlockInfo for this entry (signals it has been checked and stored!)
BlockInfo[j]=0;
}
}
if(tempdegen>1){//set of degenerate directions found so also store the reference vector and at the same time store degeneracy info!
degen.push_back(tempdegen);//store degeneracy info (note the "inverse" order of this info w.r.t order or stored vectors in the matrix)
++kcurved;//add 1 to the number of curved directions
for(size_t m=0;m!=d2V.size();++m){//Store reference degenerate direction
if(abs(MtestOut[d2V.size()*(k+kflat+1)+m])>epsilon)//Non-zero element
Mbasis[d2V.size()-(kflat+kcurved)][m]=MtestOut[d2V.size()*(k+kflat+1)+m];
else
Mbasis[d2V.size()-(kflat+kcurved)][m]=0e0; //Chop numerical errors
}
//Destroy BlockInfo for this entry (signals it has been checked and stored!)
BlockInfo[k]=0;
}
}
}
//Store non-degenerate curved eigen-directions on top
for(size_t k=0; k<d2V.size()-kflat;++k){
if(BlockInfo[k]<0){//If a curved eigen-direction...
++kcurved;
for(size_t m=0;m!=d2V.size();++m){//Store next non-degenerate direction
if(abs(MtestOut[d2V.size()*(k+kflat+1)+m])>epsilon)//Non-zero element
Mbasis[d2V.size()-(kflat+kcurved)][m]=MtestOut[d2V.size()*(k+kflat+1)+m];
else
Mbasis[d2V.size()-(kflat+kcurved)][m]=0e0; //Chop numerical errors
}
}
}
//Storing blocks now!
double kmix=0;//Loop variable over block directions...
for(size_t k=0; k<d2V.size()-kflat;++k){
if(BlockInfo[k]>0){
//Organising block;
BlockDims.push_back(1);//Initialise block dimension to 1
for(size_t m=0;m!=d2V.size();++m){//Store first direction of the block
if(abs(MtestOut[d2V.size()*(k+kflat+1)+m])>epsilon)
Mbasis[kmix][m]=MtestOut[d2V.size()*(k+kflat+1)+m];//Non-zero element
else
Mbasis[kmix][m]=0e0; //Chop numerical errors
}
++kmix;//Step to the next direction
for(size_t j=k+1; j<d2V.size()-kflat;++j){//Loop over other directions to find the ones belonging to the block
if(BlockInfo[j]==BlockInfo[k]){
++BlockDims.back();//Increse dimension of block
for(size_t m=0;m!=d2V.size();++m){//Store next direction in the block
if(abs(MtestOut[d2V.size()*(j+kflat+1)+m])>epsilon)
Mbasis[kmix][m]=MtestOut[d2V.size()*(j+kflat+1)+m];//Non-zero element
else
Mbasis[kmix][m]=0e0; //Chop numerical errors
}
++kmix;//Step to next block element
//Destroy BlockInfo for this entry (signals it has been checked and stored!)
BlockInfo[j]=0;
}
}
}
}
}
void ConvertEigen(deriv2V & d2V,deriv2V & d2Vnew,std::vector< std::vector<double> > & Mbasis,std::vector<unsigned int> & orderdV,unsigned int kflat,unsigned int kcurved,unsigned int NIC){
///Converting system to basis with a priori diagonal directions and blocks. The basis will be composed first by the "eigen-block" vectors and then the a priori eigen-vectors
vector<double> TestMag(d2V[0][0].size(),0); // Auxiliary vector to test matrix elements
///////////////////////////////////////////////////////////////
for(size_t a2=NIC;a2!=d2V[0][0].size();++a2){
//// compute the matrix elements for 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]]*Mbasis[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+=Mbasis[k][l]*Vtest[l];
}
d2Vnew[k][i][orderdV[a2]]=sum;
d2Vnew[i][k][orderdV[a2]]=sum;
TestMag[orderdV[a2]]+=abs(sum);
}
}
for(size_t i=0; i!=d2V.size(); ++i){
for(size_t k=i;k!=d2V.size();++k){
if(abs(d2Vnew[i][k][orderdV[a2]])<1e-10*TestMag[orderdV[a2]]){
d2Vnew[k][i][orderdV[a2]]=0e0;
d2Vnew[i][k][orderdV[a2]]=0e0;
}
}
}
}
unsigned int kmix=d2V.size()-kflat-kcurved;
//Chop rounding errors for eigen-directions
for(size_t a2=NIC;a2!=d2V[0][0].size();++a2){
for(size_t i=kmix; i!=d2V.size(); ++i){
if(i>=kmix+kcurved)
d2Vnew[i][i][orderdV[a2]]=0e0;// Diagonals of flat directions
for(size_t k=i+1;k!=d2V.size();++k){
d2Vnew[k][i][orderdV[a2]]=0e0; // All other
d2Vnew[i][k][orderdV[a2]]=0e0; // off diagonals
}
}
}
}
//////////////////////////////////////////
void GenMrot(double * Mdata,int dim){
//Generating random rotation matrix
RandGen r;
double *Mdata2, *Mdata3,*bdata;
Mdata2 = new double [dim*dim];// Maybe consider having this allocated from start???? and others below
Mdata3 = new double [dim*dim];
bdata = new double [dim];
for(size_t i=0;i!=dim; ++i){
for(size_t j=0;j!=dim; ++j){
//Initialise to a normal distributed random matrix
Mdata[dim*i+j]=2*r(0,1)-1e0;
if(Mdata[dim*i+j]<0)
Mdata[dim*i+j]=-sqrt(2e0*log(1e0/abs(Mdata[dim*i+j])));
else
Mdata[dim*i+j]=sqrt(2e0*log(1e0/abs(Mdata[dim*i+j])));
}
}
#ifdef DEBUG
std::clog << "Printing out Gaussian matrix used to obtain Mrot block..." << endl;
for(size_t i=0;i!=dim; ++i){
for(size_t j=0;j!=dim; ++j)
std::clog<< Mdata[dim*i+j]<< "\t";
std::clog<< endl;
}
std::clog<< endl;
#endif
gsl_matrix_view m = gsl_matrix_view_array(Mdata,dim,dim);
// exit(0);
gsl_vector_view b = gsl_vector_view_array(bdata,dim);
gsl_linalg_QR_decomp (&m.matrix, &b.vector);
gsl_matrix_view Q = gsl_matrix_view_array(Mdata2,dim,dim);
gsl_matrix_view R = gsl_matrix_view_array(Mdata3,dim,dim);
gsl_linalg_QR_unpack(&m.matrix,&b.vector, &Q.matrix, &R.matrix);
if(2*(dim/2)!=dim){
for(size_t i=0;i!=dim; ++i)
for(int j=0;j!=dim; ++j)
Mdata[dim*i+j]=gsl_matrix_get(&Q.matrix,i,j);
}
else{
for(int j=0;j!=dim; ++j){
Mdata[j]=gsl_matrix_get(&Q.matrix,1,j);
Mdata[dim+j]=gsl_matrix_get(&Q.matrix,0,j);
}
for(size_t i=2;i<dim; ++i)
for(int j=0;j!=dim; ++j)
Mdata[dim*i+j]=gsl_matrix_get(&Q.matrix,i,j);
}
delete [] Mdata2;
delete [] Mdata3;
delete [] bdata;
}
//
void ConvertWithMrot(deriv2V & d2Vnew,std::vector< std::vector<double> > & Mrot,std::vector<unsigned int> & orderdV,unsigned int kmix,unsigned int NIC){
RandGen r;
////Define a temporary matrix d2Vtemp
deriv2V d2Vtemp=d2Vnew;
//Rotate d2VnewMatrix
for(size_t a2=NIC; a2!=d2Vnew[0][0].size(); ++a2){
for(size_t i=0;i!=kmix;++i){
for(size_t j=0;j!=kmix;++j){
double sum=0e0;
for(size_t k=0;k!=kmix;++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!=kmix;++i){
for(size_t j=0;j!=kmix;++j){
double sum=0e0;
for(size_t k=0;k!=kmix;++k){
sum+=Mrot[i][k]*d2Vtemp[k][j][orderdV[a2]];
}
d2Vnew[i][j][orderdV[a2]]=sum;
}
}
}
#ifdef DEBUG
std::clog << endl<< "****** d2Vnew independent matrices in rotated basis ******" << endl;
PrintM22(d2Vnew,NIC,orderdV);
#endif
}
//////////////////////////////////////
void GenerateMrot(deriv2V & d2Vnew,std::vector< std::vector<double> > & Mrot,std::vector<unsigned int> & orderdV,unsigned int kmix,unsigned int NIC){
RandGen r;
//Generating random rotation matrix
double *Mdata, *Mdata2, *Mdata3,*bdata;
Mdata = new double [(kmix)*(kmix)];
Mdata2 = new double [(kmix)*(kmix)];
Mdata3 = new double [(kmix)*(kmix)];
bdata = new double [kmix];
for(size_t i=0;i!=kmix; ++i){
for(size_t j=0;j!=kmix; ++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[(kmix)*i+j]=Mrot[i][j];
}
}
#ifdef DEBUG
std::clog << "Printing out Gaussian matrix used to obtain Mmixing..." << endl;
PrintM1(Mrot);
#endif
size_t tempdim=kmix;
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(kmix!=0){
gsl_linalg_QR_unpack(&m.matrix,&b.vector, &Q.matrix, &R.matrix);
for(size_t i=0;i!=kmix; ++i)
for(int j=0;j!=kmix; ++j)
Mrot[i][j]=gsl_matrix_get(&Q.matrix,i,j);
if(2*(kmix/2)==kmix)
Mrot[0].swap(Mrot[1]);
#ifdef DEBUG
std::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!=kmix;++i){
for(size_t j=0;j!=kmix;++j){
double sum=0e0;
for(size_t k=0;k!=kmix;++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!=kmix;++i){
for(size_t j=0;j!=kmix;++j){
double sum=0e0;
for(size_t k=0;k!=kmix;++k){
sum+=Mrot[i][k]*d2Vtemp[k][j][orderdV[a2]];
}
d2Vnew[i][j][orderdV[a2]]=sum;
}
}
}
}
#ifdef DEBUG
std::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);
std::clog << endl << "----Eigenvalues for matrix a ="<< orderdV[a2]<< endl;
for(size_t i=0; i!=d2V.size(); ++i){
std::clog << TempOut[i]<<"\t";
}
std::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)
std::clog << TempOut[(i+1)*d2V.size()+j]<<"\t";
std::clog << endl;
}
if(abs(TempOut[d2V.size()-1])>epsilon){
size_t k;
for(size_t j=0;j!=d2V.size();++j){
if(abs(TempOut[j])>epsilon*abs(TempOut[d2V.size()-1])){
k=j;
j=d2V.size()-1;
}
}
std::clog << endl << k << " Candidate flat directions "<< endl;
}
else{
std::clog << d2V.size() << " Candidate flat directions "<< endl;
}
}
std::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 IntegratorDx(double (*Integrand)(double, void *),double * Msq)
{
gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);
double result, error;
// Assign gsl function to integrate
gsl_function F;
F.function = Integrand;
F.params = Msq;
// 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, const 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, 4:34 PM (21 h, 22 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4013428
Default Alt Text
ScannerSAux.cpp (24 KB)

Event Timeline