jiri/LA_library
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

switchable random number generators

Browse Source
This commit is contained in:
Jiri Pittner 2023-11-17 21:57:28 +01:00
parent c5a2865639
commit 578ca7bab6
10 changed files with 110 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Makefile.am
View File

@ -1,6 +1,6 @@
lib_LTLIBRARIES = libla.la lib_LTLIBRARIES = libla.la
include_HEADERS = version.h simple.h vecmat3.h quaternion.h fortran.h cuda_la.h auxstorage.h davidson.h laerror.h mat.h qsort.h vec.h bisection.h diis.h la.h noncblas.h smat.h bitvector.h fourindex.h la_traits.h nonclass.h sparsemat.h sparsesmat.h csrmat.h conjgrad.h gmres.h matexp.h permutation.h polynomial.h contfrac.h graph.h include_HEADERS = version.h simple.h vecmat3.h quaternion.h fortran.h cuda_la.h auxstorage.h davidson.h laerror.h mat.h qsort.h vec.h bisection.h diis.h la.h noncblas.h smat.h bitvector.h fourindex.h la_traits.h random.h nonclass.h sparsemat.h sparsesmat.h csrmat.h conjgrad.h gmres.h matexp.h permutation.h polynomial.h contfrac.h graph.h
libla_la_SOURCES = simple.cc quaternion.cc vecmat3.cc vec.cc mat.cc smat.cc sparsemat.cc sparsesmat.cc csrmat.cc laerror.cc noncblas.cc bitvector.cc strassen.cc nonclass.cc cuda_la.cc fourindex.cc permutation.cc polynomial.cc contfrac.cc graph.cc libla_la_SOURCES = simple.cc quaternion.cc vecmat3.cc vec.cc mat.cc smat.cc sparsemat.cc sparsesmat.cc csrmat.cc laerror.cc noncblas.cc bitvector.cc strassen.cc nonclass.cc cuda_la.cc fourindex.cc permutation.cc polynomial.cc contfrac.cc graph.cc random.cc
nodist_libla_la_SOURCES = version.cc nodist_libla_la_SOURCES = version.cc
check_PROGRAMS = t test tX check_PROGRAMS = t test tX
t_SOURCES = t.cc t2.cc t_SOURCES = t.cc t2.cc

2
la_traits.h
View File

@ -44,6 +44,7 @@
//using namespace std; //using namespace std;
#include "laerror.h" #include "laerror.h"
#include "random.h"
#include "cuda_la.h" #include "cuda_la.h"
@ -70,6 +71,7 @@ extern "C" {
namespace LA { namespace LA {
//forward declarations //forward declarations
template<typename C> class NRVec; template<typename C> class NRVec;
template<typename C> class NRVec_from1; template<typename C> class NRVec_from1;

12
mat.cc
View File

@ -1314,7 +1314,7 @@ void NRMat<double>::randomize(const double &x) {
#endif #endif
for(register int i=0; i<nn; ++i){ for(register int i=0; i<nn; ++i){
for(register int j=0; j<mm; ++j){ for(register int j=0; j<mm; ++j){
(*this)(i,j) = x*(2.*random()/(1. + RAND_MAX) - 1.); (*this)(i,j) = x*RANDDOUBLESIGNED();
} }
} }
#ifdef CUDALA #ifdef CUDALA
@ -1322,7 +1322,7 @@ void NRMat<double>::randomize(const double &x) {
NRMat<double> tmp(nn, mm, cpu); NRMat<double> tmp(nn, mm, cpu);
double *tmp_data = tmp; double *tmp_data = tmp;
for(register size_t i=0; i<(size_t)nn*mm; ++i){ for(register size_t i=0; i<(size_t)nn*mm; ++i){
tmp_data[i] = x*(2.*random()/(1. + RAND_MAX) - 1.); tmp_data[i] = x*RANDDOUBLESIGNED();
} }
tmp.moveto(this->location); tmp.moveto(this->location);
*this |= tmp; *this |= tmp;
@ -1343,8 +1343,8 @@ void NRMat<std::complex<double> >::randomize(const double &x) {
#endif #endif
for(register int i=0; i<nn; ++i){ for(register int i=0; i<nn; ++i){
for(register int j=0; j<mm; ++j){ for(register int j=0; j<mm; ++j){
const double re = x*(2.*random()/(1. + RAND_MAX) - 1.); const double re = x*RANDDOUBLESIGNED();
const double im = x*(2.*random()/(1. + RAND_MAX) - 1.); const double im = x*RANDDOUBLESIGNED();
(*this)(i,j) = std::complex<double>(re, im); (*this)(i,j) = std::complex<double>(re, im);
} }
} }
@ -1353,8 +1353,8 @@ void NRMat<std::complex<double> >::randomize(const double &x) {
NRMat<std::complex<double> > tmp(nn, mm, cpu); NRMat<std::complex<double> > tmp(nn, mm, cpu);
std::complex<double> *tmp_data = tmp; std::complex<double> *tmp_data = tmp;
for(register size_t i=0; i<(size_t)nn*mm; ++i){ for(register size_t i=0; i<(size_t)nn*mm; ++i){
const double re = x*(2.*random()/(1. + RAND_MAX) - 1.); const double re = x*RANDDOUBLESIGNED();
const double im = x*(2.*random()/(1. + RAND_MAX) - 1.); const double im = x*RANDDOUBLESIGNED();
tmp_data[i] = std::complex<double>(re, im); tmp_data[i] = std::complex<double>(re, im);
} }
tmp.moveto(this->location); tmp.moveto(this->location);

2
permutation.cc
View File

@ -167,7 +167,7 @@ this->copyonwrite();
this->identity(); this->identity();
for(int i=n-1; i>=1; --i) for(int i=n-1; i>=1; --i)
{ {
int j= random()%(i+1); int j= RANDINT32()%(i+1);
T tmp = (*this)[i+1]; T tmp = (*this)[i+1];
(*this)[i+1]=(*this)[j+1]; (*this)[i+1]=(*this)[j+1];
(*this)[j+1]=tmp; (*this)[j+1]=tmp;

8
quaternion.cc
View File

@ -417,15 +417,15 @@ void Quaternion<T>::random_rotation()
T s1,s2,s; T s1,s2,s;
do do
{ {
q[0]=2.*random()/(1. + RAND_MAX) - 1.; q[0]= RANDDOUBLESIGNED();
q[1]=2.*random()/(1. + RAND_MAX) - 1.; q[1]= RANDDOUBLESIGNED();
s1 = q[0]*q[0] + q[1]*q[1]; s1 = q[0]*q[0] + q[1]*q[1];
} }
while(s1>1); while(s1>1);
do do
{ {
q[2]=2.*random()/(1. + RAND_MAX) - 1.; q[2]= RANDDOUBLESIGNED();
q[3]=2.*random()/(1. + RAND_MAX) - 1.; q[3]= RANDDOUBLESIGNED();
s2 = q[2]*q[2] + q[3]*q[3]; s2 = q[2]*q[2] + q[3]*q[3];
} }
while(s2>1); while(s2>1);

12
random.cc Normal file
View File

@ -0,0 +1,12 @@
#include <stdlib.h>
#include "random.h"
namespace LA {
extern "C" {
WEAK_SYMBOL double randdouble() {return random()/(1.+RAND_MAX);}
WEAK_SYMBOL double randdoublesigned() {return 2.*random()/(1.+RAND_MAX)-1.;}
WEAK_SYMBOL int randint32() {return random();}
}
}//namespace

33
random.h Normal file
View File

@ -0,0 +1,33 @@
#ifndef _RANDOM_H
#define _RANDOM_H
namespace LA {
//RANDOM numbers defaulting to standard library but switchable to user's functions
#ifndef RANDDOUBLE
#define RANDDOUBLE randdouble
#endif
#ifndef RANDDOUBLESIGNED
#define RANDDOUBLESIGNED randdoublesigned
#endif
#ifndef RANDINT32
#define RANDINT32 randint32
#endif
extern "C" {
extern double randdouble();
extern double randdoublesigned();
extern int randint32();
}
#ifdef __GNUC__
#define WEAK_SYMBOL __attribute__((weak))
#else
#define WEAK_SYMBOL
#endif
}//namespace
#endif

6
smat.cc
View File

@ -249,7 +249,7 @@ void NRSMat<double>::randomize(const double &x) {
NOT_GPU(*this); NOT_GPU(*this);
for(size_t i=0; i<NN2; ++i){ for(size_t i=0; i<NN2; ++i){
v[i] = x*(2.*random()/(1.+RAND_MAX) -1.); v[i] = x*RANDDOUBLESIGNED();
} }
} }
@ -259,8 +259,8 @@ void NRSMat<double>::randomize(const double &x) {
******************************************************************************/ ******************************************************************************/
template<> template<>
void NRSMat<std::complex<double> >::randomize(const double &x) { void NRSMat<std::complex<double> >::randomize(const double &x) {
for(register size_t i=0; i<NN2; ++i) v[i].real(x*(2.*random()/(1. + RAND_MAX) -1.)); for(register size_t i=0; i<NN2; ++i) v[i].real(x*RANDDOUBLESIGNED());
for(register size_t i=0; i<NN2; ++i) v[i].imag(x*(2.*random()/(1. + RAND_MAX) -1.)); for(register size_t i=0; i<NN2; ++i) v[i].imag(x*RANDDOUBLESIGNED());
for(register int i=0; i<nn; ++i){ for(register int i=0; i<nn; ++i){
for(register int j=0; j<=i; ++j){ for(register int j=0; j<=i; ++j){
if(i == j) v[i*(size_t)(i+1)/2+j].imag(0.); //hermitean if(i == j) v[i*(size_t)(i+1)/2+j].imag(0.); //hermitean

90
t.cc
View File

@ -54,7 +54,7 @@ void f2(double *c)
inline int randind(const int n) inline int randind(const int n)
{ {
return int(random()/(1.+RAND_MAX)*n); return int(RANDDOUBLE()*n);
} }
complex<double> mycident (const complex<double>&x) {return x;} complex<double> mycident (const complex<double>&x) {return x;}
@ -640,7 +640,7 @@ for(int n=8; n<=1024*1024;n+=n)
cout << "\n\n\ntiming for size "<<n<<endl; cout << "\n\n\ntiming for size "<<n<<endl;
if(n<=512) { if(n<=512) {
NRMat<double> a(0.,n,n); NRMat<double> a(0.,n,n);
for(int i=0; i<sparsity;i++) a(randind(n),randind(n))=random()/(1.+RAND_MAX); for(int i=0; i<sparsity;i++) a(randind(n),randind(n))=RANDDOUBLE();
double t0=clock()/((double) (CLOCKS_PER_SEC)); double t0=clock()/((double) (CLOCKS_PER_SEC));
//cout <<a; //cout <<a;
NRMat<double> b(exp(a)); NRMat<double> b(exp(a));
@ -653,7 +653,7 @@ for(int n=8; n<=1024*1024;n+=n)
} }
else else
{ {
for(int i=0; i<sparsity;i++) aa.add(randind(n),randind(n),random()/(1.+RAND_MAX)); for(int i=0; i<sparsity;i++) aa.add(randind(n),randind(n),RANDDOUBLE());
} }
//cout <<aa; //cout <<aa;
double t2=clock()/((double) (CLOCKS_PER_SEC)); double t2=clock()/((double) (CLOCKS_PER_SEC));
@ -672,10 +672,10 @@ if(0)
int n; int n;
cin>>n; cin>>n;
SparseMat<double> aa(n,n); SparseMat<double> aa(n,n);
for(int i=0; i<sparsity;i++) aa.add(randind(n),randind(n),random()/(1.+RAND_MAX)); for(int i=0; i<sparsity;i++) aa.add(randind(n),randind(n),RANDDOUBLE());
SparseMat<double> bb=exp(aa); SparseMat<double> bb=exp(aa);
NRVec<double> v(n); NRVec<double> v(n);
for(int i=0; i<n;++i) v[i]=random()/(1.+RAND_MAX); for(int i=0; i<n;++i) v[i]=RANDDOUBLE();
NRVec<double> res1=bb*v; NRVec<double> res1=bb*v;
NRVec<double> res2=exptimes(aa,v); NRVec<double> res2=exptimes(aa,v);
cout <<"difference = "<<(res1-res2).norm()<<endl; cout <<"difference = "<<(res1-res2).norm()<<endl;
@ -687,8 +687,8 @@ int n,k,m;
cin >>n>>k>>m; cin >>n>>k>>m;
{ {
NRMat<double> a(n,k),b(k,m),c(n,m),d(n,m); NRMat<double> a(n,k),b(k,m),c(n,m),d(n,m);
for(int i=0;i<n;++i) for(int j=0;j<k;++j) a(i,j)= random()/(1.+RAND_MAX); for(int i=0;i<n;++i) for(int j=0;j<k;++j) a(i,j)= RANDDOUBLE();
for(int i=0;i<k;++i) for(int j=0;j<m;++j) b(i,j)= random()/(1.+RAND_MAX); for(int i=0;i<k;++i) for(int j=0;j<m;++j) b(i,j)= RANDDOUBLE();
c.gemm(0., a, 'n', b, 'n', .6); c.gemm(0., a, 'n', b, 'n', .6);
SparseMat<double> aa(a); SparseMat<double> aa(a);
d.gemm(0., aa, 'n', b, 'n', .6); d.gemm(0., aa, 'n', b, 'n', .6);
@ -697,8 +697,8 @@ cout <<"test error = "<<(c-d).norm()<<endl;
} }
{ {
NRMat<double> a(k,n),b(k,m),c(n,m),d(n,m); NRMat<double> a(k,n),b(k,m),c(n,m),d(n,m);
for(int i=0;i<k;++i) for(int j=0;j<n;++j) a(i,j)= random()/(1.+RAND_MAX); for(int i=0;i<k;++i) for(int j=0;j<n;++j) a(i,j)= RANDDOUBLE();
for(int i=0;i<k;++i) for(int j=0;j<m;++j) b(i,j)= random()/(1.+RAND_MAX); for(int i=0;i<k;++i) for(int j=0;j<m;++j) b(i,j)= RANDDOUBLE();
c.gemm(0., a, 't', b, 'n', .7); c.gemm(0., a, 't', b, 'n', .7);
SparseMat<double> aa(a); SparseMat<double> aa(a);
d.gemm(0., aa, 't', b, 'n', .7); d.gemm(0., aa, 't', b, 'n', .7);
@ -707,8 +707,8 @@ cout <<"test error = "<<(c-d).norm()<<endl;
} }
{ {
NRMat<double> a(n,k),b(m,k),c(n,m),d(n,m); NRMat<double> a(n,k),b(m,k),c(n,m),d(n,m);
for(int i=0;i<n;++i) for(int j=0;j<k;++j) a(i,j)= random()/(1.+RAND_MAX); for(int i=0;i<n;++i) for(int j=0;j<k;++j) a(i,j)= RANDDOUBLE();
for(int i=0;i<m;++i) for(int j=0;j<k;++j) b(i,j)= random()/(1.+RAND_MAX); for(int i=0;i<m;++i) for(int j=0;j<k;++j) b(i,j)= RANDDOUBLE();
c.gemm(0., a, 'n', b, 't', .8); c.gemm(0., a, 'n', b, 't', .8);
SparseMat<double> aa(a); SparseMat<double> aa(a);
d.gemm(0., aa, 'n', b, 't', .8); d.gemm(0., aa, 'n', b, 't', .8);
@ -717,8 +717,8 @@ cout <<"test error = "<<(c-d).norm()<<endl;
} }
{ {
NRMat<double> a(k,n),b(m,k),c(n,m),d(n,m); NRMat<double> a(k,n),b(m,k),c(n,m),d(n,m);
for(int i=0;i<k;++i) for(int j=0;j<n;++j) a(i,j)= random()/(1.+RAND_MAX); for(int i=0;i<k;++i) for(int j=0;j<n;++j) a(i,j)= RANDDOUBLE();
for(int i=0;i<m;++i) for(int j=0;j<k;++j) b(i,j)= random()/(1.+RAND_MAX); for(int i=0;i<m;++i) for(int j=0;j<k;++j) b(i,j)= RANDDOUBLE();
c.gemm(0., a, 't', b, 't', .9); c.gemm(0., a, 't', b, 't', .9);
SparseMat<double> aa(a); SparseMat<double> aa(a);
d.gemm(0., aa, 't', b, 't', .9); d.gemm(0., aa, 't', b, 't', .9);
@ -764,8 +764,8 @@ cin >>n;
NRMat<double> a(n,n); NRMat<double> a(n,n);
for(int i=0;i<n;++i) for(int j=0;j<=i;++j) for(int i=0;i<n;++i) for(int j=0;j<=i;++j)
{ {
a(i,j)= random()/(1.+RAND_MAX); a(i,j)= RANDDOUBLE();
a(j,i)= random()/(1.+RAND_MAX); a(j,i)= RANDDOUBLE();
} }
NRMat<double> b; b|=a; NRMat<double> b; b|=a;
NRVec<double> er(n),ei(n); NRVec<double> er(n),ei(n);
@ -810,7 +810,7 @@ cin>>n;
NRMat<double> a(n,n); NRMat<double> a(n,n);
for(int i=0;i<n;++i) for(int j=0;j<i;++j) for(int i=0;i<n;++i) for(int j=0;j<i;++j)
{ {
a(i,j)= random()/(1.+RAND_MAX); a(i,j)= RANDDOUBLE();
a(j,i)= -a(i,j); a(j,i)= -a(i,j);
} }
cout <<"a matrix \n"<<a; cout <<"a matrix \n"<<a;
@ -836,9 +836,9 @@ NRMat<double> a(n,n);
NRVec<double>u(n),v,w; NRVec<double>u(n),v,w;
for(int i=0;i<n;++i) for(int i=0;i<n;++i)
{ {
u[i]=f*random()/(1.+RAND_MAX); u[i]=f*RANDDOUBLE();
for(int j=0;j<n;++j) for(int j=0;j<n;++j)
a(i,j)= f*random()/(1.+RAND_MAX); a(i,j)= f*RANDDOUBLE();
} }
//cout <<"a matrix \n"<<a; //cout <<"a matrix \n"<<a;
//cout<<"EXP\n"; //cout<<"EXP\n";
@ -873,7 +873,7 @@ cin>>n;
NRMat<double> a(n,n); NRMat<double> a(n,n);
for(int i=0;i<n;++i) for(int j=0;j<=i;++j) for(int i=0;i<n;++i) for(int j=0;j<=i;++j)
{ {
a(i,j)= .1*random()/(1.+RAND_MAX); a(i,j)= .1*RANDDOUBLE();
a(j,i)= a(i,j); a(j,i)= a(i,j);
} }
NRMat<double> b=exp(a); NRMat<double> b=exp(a);
@ -1029,7 +1029,7 @@ cin >>n;
SparseMat<double> a(n,n); SparseMat<double> a(n,n);
for(int i=0;i<n;++i) for(int j=0;j<=i;++j) for(int i=0;i<n;++i) for(int j=0;j<=i;++j)
{ {
a.add(i,j,random()/(1.+RAND_MAX)); a.add(i,j,RANDDOUBLE());
} }
a.setsymmetric(); a.setsymmetric();
NRSMat<double> aa(a); NRSMat<double> aa(a);
@ -1061,11 +1061,11 @@ cin >>n;
NRSMat<double> a(n); NRSMat<double> a(n);
for(int i=0;i<n;++i) for(int j=0;j<=i;++j) for(int i=0;i<n;++i) for(int j=0;j<=i;++j)
{ {
a(j,i)=a(i,j)=random()/(1.+RAND_MAX)*(i==j?10.:1.); a(j,i)=a(i,j)=RANDDOUBLE()*(i==j?10.:1.);
} }
cout <<a; cout <<a;
NRMat<double> y(1,n); NRMat<double> y(1,n);
for(int i=0;i<n;++i) y(0,i)=random()/(1.+RAND_MAX); for(int i=0;i<n;++i) y(0,i)=RANDDOUBLE();
cout <<y; cout <<y;
linear_solve(a,&y,&d); linear_solve(a,&y,&d);
cout << y; cout << y;
@ -1078,7 +1078,7 @@ int n;
cin >>n; cin >>n;
SparseMat<double> a(n,n); SparseMat<double> a(n,n);
int spars=n*n/3; int spars=n*n/3;
for(int i=0; i<spars;i++) a.add(randind(n),randind(n),random()/(1.+RAND_MAX)); for(int i=0; i<spars;i++) a.add(randind(n),randind(n),RANDDOUBLE());
NRMat<double> aa(a); NRMat<double> aa(a);
NRVec<double> v(aa[0],n*n); NRVec<double> v(aa[0],n*n);
@ -1115,7 +1115,7 @@ NRVec<double> rr(n);
for(int i=0;i<n;++i) for(int j=0;j<=i;++j) for(int i=0;i<n;++i) for(int j=0;j<=i;++j)
{ {
a(i,j)= random()/(1.+RAND_MAX); a(i,j)= RANDDOUBLE();
if(i==j) a(i,i)+= .5*(i-n); if(i==j) a(i,i)+= .5*(i-n);
} }
@ -1151,8 +1151,8 @@ NRVec<double> rr(n),ii(n);
double tmp=0.; double tmp=0.;
for(int i=0;i<n;++i) for(int j=0;j<n;++j) for(int i=0;i<n;++i) for(int j=0;j<n;++j)
{ {
a(i,j)= random()/(1.+RAND_MAX); a(i,j)= RANDDOUBLE();
a(j,i)= random()/(1.+RAND_MAX); a(j,i)= RANDDOUBLE();
if(i==j) a(i,i)+= .5*(i-n); if(i==j) a(i,i)+= .5*(i-n);
tmp+= (a(i,j)-a(j,i))*(a(i,j)-a(j,i)); tmp+= (a(i,j)-a(j,i))*(a(i,j)-a(j,i));
} }
@ -1190,8 +1190,8 @@ int n,m;
cin >>n>>m; cin >>n>>m;
SparseMat<double> aa(n,n); SparseMat<double> aa(n,n);
aa.setsymmetric(); aa.setsymmetric();
for(int i=0; i<sparsity;i++) aa.add(randind(n),randind(n),random()/(1.+RAND_MAX)); for(int i=0; i<sparsity;i++) aa.add(randind(n),randind(n),RANDDOUBLE());
for(int i=0; i<n; ++i) aa.add(i,i,500*random()/(1.+RAND_MAX)); for(int i=0; i<n; ++i) aa.add(i,i,500*RANDDOUBLE());
NRVec<double> r(m); NRVec<double> r(m);
davidson(aa,r,(NRVec<double> *)NULL,"eivecs",m,1,1e-5,0,300,300); davidson(aa,r,(NRVec<double> *)NULL,"eivecs",m,1,1e-5,0,300,300);
cout <<r; cout <<r;
@ -1207,13 +1207,13 @@ int n,m;
cin >>n>>m; cin >>n>>m;
SparseMat<double> aa(n,n); SparseMat<double> aa(n,n);
aa.setsymmetric(); aa.setsymmetric();
for(int i=0; i<sparsity;i++) aa.add(randind(n),randind(n),random()/(1.+RAND_MAX)); for(int i=0; i<sparsity;i++) aa.add(randind(n),randind(n),RANDDOUBLE());
for(int i=0; i<n; ++i) aa.add(i,i,500*random()/(1.+RAND_MAX)); for(int i=0; i<n; ++i) aa.add(i,i,500*RANDDOUBLE());
NRVec<double> r(m); NRVec<double> r(m);
NRVec<double> r2(m); NRVec<double> r2(m);
davidson(aa,r,(NRVec<double> *)NULL,"eivecs",m,1,1e-5,1,300,300); davidson(aa,r,(NRVec<double> *)NULL,"eivecs",m,1,1e-5,1,300,300);
SparseMat<double> bb(n,n); SparseMat<double> bb(n,n);
for(int i=0; i<sparsity2;i++) bb.add(randind(n),randind(n),random()/(1.+RAND_MAX)); for(int i=0; i<sparsity2;i++) bb.add(randind(n),randind(n),RANDDOUBLE());
SparseMat<double> e1,e2,cc; SparseMat<double> e1,e2,cc;
e1=exp(bb); e1=exp(bb);
e2=exp(bb*-1.); e2=exp(bb*-1.);
@ -1236,12 +1236,12 @@ cin >>n>>m;
{ {
int k= randind(n); int k= randind(n);
int l= randind(n); int l= randind(n);
double a=random()/(1.+RAND_MAX); double a=RANDDOUBLE();
double b=random()/(1.+RAND_MAX)-.5; double b=RANDDOUBLE()-.5;
aa.add(k,l,a); aa.add(k,l,a);
aa.add(l,k,a+b/20); aa.add(l,k,a+b/20);
} }
for(int i=0; i<n; ++i) aa.add(i,i,500*random()/(1.+RAND_MAX)); for(int i=0; i<n; ++i) aa.add(i,i,500*RANDDOUBLE());
NRVec<double> r(m); NRVec<double> r(m);
davidson(aa,r,(NRVec<double> *)NULL,"eivecs",m,1,1e-5,0,300,300); davidson(aa,r,(NRVec<double> *)NULL,"eivecs",m,1,1e-5,0,300,300);
cout <<r; cout <<r;
@ -1259,7 +1259,7 @@ NRVec<double> x(m);
for(int i=0;i<n;++i) for(int j=0;j<m;++j) for(int i=0;i<n;++i) for(int j=0;j<m;++j)
{ {
//a(j,i)= 2*i*i*i-5*j+ +9*8*7*6*5*4*3*2/(i+j+1.)+3*(i*i+2*j*j*j); //a(j,i)= 2*i*i*i-5*j+ +9*8*7*6*5*4*3*2/(i+j+1.)+3*(i*i+2*j*j*j);
a(i,j)= random()/(1.+RAND_MAX); a(i,j)= RANDDOUBLE();
if(i==j) a(i,i)+= .5*(i-n); if(i==j) a(i,i)+= .5*(i-n);
} }
for(int i=0;i<n;++i) b[i] = i; for(int i=0;i<n;++i) b[i] = i;
@ -1297,9 +1297,9 @@ NRVec<double> b(n);
NRVec<double> x(m); NRVec<double> x(m);
//tridiagonal //tridiagonal
for(int i=0; i<n; ++i) aa.add(i,i,random()/(1.+RAND_MAX)); for(int i=0; i<n; ++i) aa.add(i,i,RANDDOUBLE());
for(int i=0; i<n-1; i+=1) aa.add(i,i+1,.002*random()/(1.+RAND_MAX)); for(int i=0; i<n-1; i+=1) aa.add(i,i+1,.002*RANDDOUBLE());
for(int i=0; i<n-1; i+=1) aa.add(i+1,i,.002*random()/(1.+RAND_MAX)); for(int i=0; i<n-1; i+=1) aa.add(i+1,i,.002*RANDDOUBLE());
for(int i=0;i<n;++i) b[i] = i+1; for(int i=0;i<n;++i) b[i] = i+1;
gmres(aa,b,x,1,1e-20,20,1,1,1,1000,1); gmres(aa,b,x,1,1e-20,20,1,1,1,1000,1);
@ -1323,7 +1323,7 @@ int dim;
cin>>dim; cin>>dim;
NRVec<double> solution(dim), deviation(dim); NRVec<double> solution(dim), deviation(dim);
for(i=0; i<dim; ++i) solution[i]=i&1 ? i/2.:-i-3.; for(i=0; i<dim; ++i) solution[i]=i&1 ? i/2.:-i-3.;
for(i=0; i<dim; ++i) deviation[i]= (i&2 ? 1:-1) * random()/(1.+RAND_MAX); for(i=0; i<dim; ++i) deviation[i]= (i&2 ? 1:-1) * RANDDOUBLE();
double norm=1e100; double norm=1e100;
for(int iter=1; iter<100 && norm>1e-8 ; ++iter) for(int iter=1; iter<100 && norm>1e-8 ; ++iter)
{ {
@ -1400,7 +1400,7 @@ cout <<"test eigenvalues "<<hlp.norm()<<endl;
NRMat<double> ader(n,n); NRMat<double> ader(n,n);
for(int i=0; i<n; ++i) for(int i=0; i<n; ++i)
for(int j=0; j<n; ++j) for(int j=0; j<n; ++j)
ader(i,j) = 2.*random()/(1.+RAND_MAX) -1.; ader(i,j) = RANDDOUBLESIGNED();
cout <<"eigenvalues\n"<<wr<<endl; cout <<"eigenvalues\n"<<wr<<endl;
@ -1722,7 +1722,7 @@ int n;
cin >>n; cin >>n;
NRMat<double> bb(0.,n,n); NRMat<double> bb(0.,n,n);
int nn=0; int nn=0;
for(int i=0; i<n; ++i) for(int j=i; j<n; ++j) {if((double)random()/RAND_MAX>0.995 || i==j) {++nn; bb(i,j)=bb(j,i)=(double)random()/RAND_MAX*(i==j?10.:.5/(i+j)*(random()>RAND_MAX/2?1:-1));}} for(int i=0; i<n; ++i) for(int j=i; j<n; ++j) {if(RANDDOUBLE()>0.995 || i==j) {++nn; bb(i,j)=bb(j,i)=RANDDOUBLE()*(i==j?10.:.5/(i+j)*(RANDDOUBLE()>.5?1:-1));}}
bb=bb*bb.transpose(); bb=bb*bb.transpose();
//cout <<bb; //cout <<bb;
nn=0; nn=0;
@ -1786,8 +1786,8 @@ if(0)
int n; int n;
cin >>n; cin >>n;
SparseSMat<double> bh(n,n); SparseSMat<double> bh(n,n);
for(int i=0; i<=n/400; ++i) for(int j=i; j<n; ++j) {if((double)random()/RAND_MAX>0.995 || i==j) for(int i=0; i<=n/400; ++i) for(int j=i; j<n; ++j) {if(RANDDOUBLE()>0.995 || i==j)
{bh.add(i,j,(double)random()/RAND_MAX*(i==j?10.:(random()>RAND_MAX/2?1:-1)),false);}} {bh.add(i,j,RANDDOUBLE()*(i==j?10.:(RANDDOUBLE()>.5?1:-1)),false);}}
if(n<1000) cout <<"Random matrix\n"<<bh; if(n<1000) cout <<"Random matrix\n"<<bh;
SparseSMat<double> bb(n,n); SparseSMat<double> bb(n,n);
bb.gemm(0.,bh,'c',bh,'n',1.); bb.gemm(0.,bh,'c',bh,'n',1.);
@ -2440,7 +2440,7 @@ simple_linfit<double,3> fit;
double funcs[3]; double funcs[3];
for(int i=0; i<100; ++i) for(int i=0; i<100; ++i)
{ {
double x = 10*(2.*random()/(1. + RAND_MAX) - 1.); double x = 10*RANDDOUBLESIGNED();
funcs[0]=1; funcs[0]=1;
funcs[1]=x; funcs[1]=x;
funcs[2]=x*x; funcs[2]=x*x;
@ -2585,7 +2585,7 @@ if(sizeof(int)!=read(f,&seed,sizeof(int))) laerror("cannot read /dev/random");
close(f); close(f);
srand(seed); srand(seed);
double t=random(); double t=RANDINT32();
t/=(1L<<30); t/=(1L<<30);
t-=1; t-=1;
cout <<"theta = "<<t<<endl; cout <<"theta = "<<t<<endl;

4
vec.cc
View File

@ -194,7 +194,7 @@ void NRVec<double>::randomize(const double &x){
NOT_GPU(*this); NOT_GPU(*this);
for(register int i=0; i<nn; ++i){ for(register int i=0; i<nn; ++i){
v[i] = x*(2.*random()/(1. + RAND_MAX) - 1.); v[i] = x*RANDDOUBLESIGNED();
} }
} }
@ -211,7 +211,7 @@ void NRVec<std::complex<double> >::randomize(const double &x) {
NOT_GPU(*this); NOT_GPU(*this);
for(register int i=0; i<nn; ++i){ for(register int i=0; i<nn; ++i){
v[i] = std::complex<double>(x*(2.*random()/(1. + RAND_MAX) - 1.), x*(2.*random()/(1. + RAND_MAX) - 1.)); v[i] = std::complex<double>(x*RANDDOUBLESIGNED(), x*RANDDOUBLESIGNED());
} }
} }