CoreOther.cpp

Go to the documentation of this file.

//Copyright (c) 2004-2005, Baris Sumengen
//All rights reserved.
//
// CIMPL Matrix Performance Library
//
//Redistribution and use in source and binary
//forms, with or without modification, are
//permitted provided that the following
//conditions are met:
//
//    * No commercial use is allowed. 
//    This software can only be used
//    for non-commercial purposes. This 
//    distribution is mainly intended for
//    academic research and teaching.
//    * Redistributions of source code must
//    retain the above copyright notice, this
//    list of conditions and the following
//    disclaimer.
//    * Redistributions of binary form must
//    mention the above copyright notice, this
//    list of conditions and the following
//    disclaimer in a clearly visible part 
//    in associated product manual, 
//    readme, and web site of the redistributed 
//    software.
//    * Redistributions in binary form must
//    reproduce the above copyright notice,
//    this list of conditions and the
//    following disclaimer in the
//    documentation and/or other materials
//    provided with the distribution.
//    * The name of Baris Sumengen may not be
//    used to endorse or promote products
//    derived from this software without
//    specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
//HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
//EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
//NOT LIMITED TO, THE IMPLIED WARRANTIES OF
//MERCHANTABILITY AND FITNESS FOR A PARTICULAR
//PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
//CONTRIBUTORS BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
//EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
//OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
//DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
//HOWEVER CAUSED AND ON ANY THEORY OF
//LIABILITY, WHETHER IN CONTRACT, STRICT
//LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
//OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
//OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//POSSIBILITY OF SUCH DAMAGE.



#include "./Core.h"


namespace MathCore
{
        void Ind2Sub(int rows, int cols, Vector<int>& ind, Vector<int>& I, Vector<int>& J)
        {
                I = Vector<int>(ind.Length());
                J = Vector<int>(ind.Length());
                
                int maxElem = rows*cols;

                for(int i=0; i<ind.Length(); i++)
                {
                        int elem = ind.ElemNC(i);
                        if(elem>=maxElem)
                        {
                                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                                Utility::RunTimeError("Index is larger than the size of the matrix!");
                        }
                        I.ElemNC(i) = elem % rows;
                        J.ElemNC(i) = elem / rows;
                }
        }
        
        Vector<int> Sub2Ind(int rows, int cols, Vector<int>& I, Vector<int>& J)
        {
                if(I.Length() != J.Length())
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Index lengths are not the same!");
                }
                
                Vector<int> temp(I.Length());
                
                for(int i=0; i<I.Length(); i++)
                {
                        if(I.ElemNC(i)>=rows || J.ElemNC(i)>=cols)
                        {
                                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                                Utility::RunTimeError("Index is larger than the dimensions of the matrix!");
                        }
                        temp.ElemNC(i) = I.ElemNC(i) + rows*J.ElemNC(i);
                }
                
                return temp;
        }


        
        Vector<double> LinSpace(double x, double y, int N)
        {
                double step = (y-x)/(N-1);
                Vector<double> temp(N);
                temp.ElemNC(0) = x;
                for(int i=1; i<N-1; i++)
                {
                        temp.ElemNC(i) = temp.ElemNC(i-1) + step;
                }
                temp.ElemNC(N-1) = y;
                return temp;
        }
        
        Vector<double> LinSpace(double x, double y)
        {
                return LinSpace(x,y,100);
        }


        Vector<double> LogSpace(double x, double y, int N)
        {
                double step = (y-x)/(N-1);
                Vector<double> temp(N);
                temp.ElemNC(0) = pow(10.0,x);
                for(int i=1; i<N-1; i++)
                {
                        temp.ElemNC(i) = pow(10.0, x+i*step);
                }
                temp.ElemNC(N-1) = pow(10.0,y);
                return temp;
        }
        Vector<double> LogSpace(double x, double y)
        {
                return LogSpace(x,y,50);
        }

        
        Matrix<float>& RandN(Matrix<float>& m)
        {
                if(!RandomGen::Initialized())
                {
                        RandomGen::Initialize();
                }
                
                int le = m.Length()/2;
                int rem = m.Length() % 2;
                float rsq;
                float v1;
                float v2;
                for(int i=0;i<le;i++)
                {
                        do 
                        {
                                float r1 = rand()/(float)RAND_MAX;
                                v1 = 2.0f*r1 - 1.0f;
                                float r2 = rand()/(float)RAND_MAX;
                                v2 = 2.0f*r2 - 1.0f;
                                rsq=v1*v1+v2*v2;
                        } 
                        while (rsq >= 1.0 || rsq == 0.0);
                        float fac = (float)sqrt(-2.0*log(rsq)/rsq);
                        m.ElemNC(2*i) = v1*fac;
                        m.ElemNC(2*i+1) = v2*fac;
                }
                if(rem == 1)
                {
                        do 
                        {
                                float r1 = rand()/(float)RAND_MAX;
                                v1 = 2.0f*r1 - 1.0f;
                                float r2 = rand()/(float)RAND_MAX;
                                v2 = 2.0f*r2 - 1.0f;
                                rsq=v1*v1+v2*v2;
                        } 
                        while (rsq >= 1.0 || rsq == 0.0);
                        float fac = (float)sqrt(-2.0*log(rsq)/rsq);
                        m.ElemNC(m.Length()-1) = v1*fac;
                }
                
                return m;
        }
        Vector<float>& RandN(Vector<float>& m)
        {
                
                if(!RandomGen::Initialized())
                {
                        RandomGen::Initialize();
                }
                
                int le = m.Length()/2;
                int rem = m.Length() % 2;
                float rsq;
                float v1;
                float v2;
                for(int i=0;i<le;i++)
                {
                        do 
                        {
                                float r1 = rand()/(float)RAND_MAX;
                                v1 = 2.0f*r1 - 1.0f;
                                float r2 = rand()/(float)RAND_MAX;
                                v2 = 2.0f*r2 - 1.0f;
                                rsq=v1*v1+v2*v2;
                        } 
                        while (rsq >= 1.0 || rsq == 0.0);
                        float fac = (float)sqrt(-2.0*log(rsq)/rsq);
                        m.ElemNC(2*i) = v1*fac;
                        m.ElemNC(2*i+1) = v2*fac;
                }
                if(rem == 1)
                {
                        do 
                        {
                                float r1 = rand()/(float)RAND_MAX;
                                v1 = 2.0f*r1 - 1.0f;
                                float r2 = rand()/(float)RAND_MAX;
                                v2 = 2.0f*r2 - 1.0f;
                                rsq=v1*v1+v2*v2;
                        } 
                        while (rsq >= 1.0 || rsq == 0.0);
                        float fac = (float)sqrt(-2.0*log(rsq)/rsq);
                        m.ElemNC(m.Length()-1) = v1*fac;
                }
                return m;
        }


        Matrix<double>& RandN(Matrix<double>& m)
        {
                if(!RandomGen::Initialized())
                {
                        RandomGen::Initialize();
                }
                
                int le = m.Length()/2;
                int rem = m.Length() % 2;
                double rsq;
                double v1;
                double v2;
                for(int i=0;i<le;i++)
                {
                        do 
                        {
                                double r1 = rand()/(double)RAND_MAX;
                                v1 = 2.0*r1 - 1.0;
                                double r2 = rand()/(double)RAND_MAX;
                                v2 = 2.0*r2 - 1.0;
                                rsq=v1*v1+v2*v2;
                        } 
                        while (rsq >= 1.0 || rsq == 0.0);
                        double fac=sqrt(-2.0*log(rsq)/rsq);
                        m.ElemNC(2*i) = v1*fac;
                        m.ElemNC(2*i+1) = v2*fac;
                }
                if(rem == 1)
                {
                        do 
                        {
                                double r1 = rand()/(double)RAND_MAX;
                                v1 = 2.0*r1 - 1.0;
                                double r2 = rand()/(double)RAND_MAX;
                                v2 = 2.0*r2 - 1.0;
                                rsq=v1*v1+v2*v2;
                        } 
                        while (rsq >= 1.0 || rsq == 0.0);
                        double fac=sqrt(-2.0*log(rsq)/rsq);
                        m.ElemNC(m.Length()-1) = v1*fac;
                }
                
                return m;
        }
        
        
        Vector<double>& RandN(Vector<double>& m)
        {
                
                if(!RandomGen::Initialized())
                {
                        RandomGen::Initialize();
                }
                
                int le = m.Length()/2;
                int rem = m.Length() % 2;
                double rsq;
                double v1;
                double v2;
                for(int i=0;i<le;i++)
                {
                        do 
                        {
                                double r1 = rand()/(double)RAND_MAX;
                                v1 = 2.0*r1 - 1.0;
                                double r2 = rand()/(double)RAND_MAX;
                                v2 = 2.0*r2 - 1.0;
                                rsq=v1*v1+v2*v2;
                        } 
                        while (rsq >= 1.0 || rsq == 0.0);
                        double fac=sqrt(-2.0*log(rsq)/rsq);
                        m.ElemNC(2*i) = v1*fac;
                        m.ElemNC(2*i+1) = v2*fac;
                }
                if(rem == 1)
                {
                        do 
                        {
                                double r1 = rand()/(double)RAND_MAX;
                                v1 = 2.0*r1 - 1.0;
                                double r2 = rand()/(double)RAND_MAX;
                                v2 = 2.0*r2 - 1.0;
                                rsq=v1*v1+v2*v2;
                        } 
                        while (rsq >= 1.0 || rsq == 0.0);
                        double fac=sqrt(-2.0*log(rsq)/rsq);
                        m.ElemNC(m.Length()-1) = v1*fac;
                }
                return m;
        }



        
        int NextPow2(int n)
        {
                double l = log((double)n)/log(2.0);
                return (int)pow(2,(int)ceil(l));
        }




        Matrix<ComplexFloat>& Rand(Matrix<ComplexFloat>& m, float max)
        {
                if(!RandomGen::Initialized())
                {
                        RandomGen::Initialize();
                }
                for(int i=0;i<m.Length();i++)
                {
                        m.ElemNC(i) = ComplexFloat( (float)(rand()/(double)RAND_MAX*max), (float)(rand()/(double)RAND_MAX*max) );
                }
                return m;
        }
        
        Vector<ComplexFloat>& Rand(Vector<ComplexFloat>& m, float max)
        {
                if(!RandomGen::Initialized())
                {
                        RandomGen::Initialize();
                }
                for(int i=0;i<m.Length();i++)
                {
                        m.ElemNC(i) = ComplexFloat( (float)(rand()/(double)RAND_MAX*max), (float)(rand()/(double)RAND_MAX*max) );
                }
                return m;
        }


        Matrix<ComplexDouble>& Rand(Matrix<ComplexDouble>& m, double max)
        {
                if(!RandomGen::Initialized())
                {
                        RandomGen::Initialize();
                }
                for(int i=0;i<m.Length();i++)
                {
                        m.ElemNC(i) = ComplexDouble( (double)(rand()/(double)RAND_MAX*max), (double)(rand()/(double)RAND_MAX*max) );
                }
                return m;
        }

        Vector<ComplexDouble>& Rand(Vector<ComplexDouble>& m, double max)
        {
                if(!RandomGen::Initialized())
                {
                        RandomGen::Initialize();
                }
                for(int i=0;i<m.Length();i++)
                {
                        m.ElemNC(i) = ComplexDouble( (double)(rand()/(double)RAND_MAX*max), (double)(rand()/(double)RAND_MAX*max) );
                }
                return m;
        }












};

---