Changeset 364ccf in git for Singular/svd_si.h


Ignore:
Timestamp:
Sep 21, 2017, 3:33:15 PM (7 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
f166afac68e9933694b0ff29efc419fc621d35c4
Parents:
2206a04585b4cd3234b11b58891ba713a1955554
Message:
fix: throw/catch in svd_si.h
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/svd_si.h

    r2206a04 r364ccf  
    3434
    3535/********************************************************************
    36 This symbol is used for debugging. Do not define it and do not remove 
    37 comments. 
     36This symbol is used for debugging. Do not define it and do not remove
     37comments.
    3838********************************************************************/
    3939//#define UNSAFE_MEM_COPY
    … …  
    130130
    131131/********************************************************************
    132 Template defining vector in memory. It is used by the basic 
     132Template defining vector in memory. It is used by the basic
    133133subroutines of linear algebra.
    134134
    135 Vector consists of Length elements of type T, starting from an element, 
    136 which Data is pointed to. Interval between adjacent elements equals 
     135Vector consists of Length elements of type T, starting from an element,
     136which Data is pointed to. Interval between adjacent elements equals
    137137the value of Step.
    138138
    … …  
    164164It is used by the basic subroutines of linear algebra.
    165165
    166 Vector consists of Length elements of type T, starting from an element, 
    167 which Data is pointed to. Interval between adjacent elements equals 
     166Vector consists of Length elements of type T, starting from an element,
     167which Data is pointed to. Interval between adjacent elements equals
    168168the value of Step.
    169169
    … …  
    694694    };
    695695
    696    
     696
    697697    const template_1d_array& operator=(const template_1d_array &rhs)
    698698    {
    … …  
    720720    };
    721721
    722    
     722
    723723    const T& operator()(int i) const
    724724    {
    … …  
    729729    };
    730730
    731    
     731
    732732    T& operator()(int i)
    733733    {
    … …  
    738738    };
    739739
    740    
     740
    741741    void setbounds( int iLow, int iHigh )
    742742    {
    … …  
    749749    };
    750750
    751    
     751
    752752    void setcontent( int iLow, int iHigh, const T *pContent )
    753753    {
    … …  
    757757    };
    758758
    759    
     759
    760760    T* getcontent()
    761761    {
    … …  
    768768    };
    769769
    770    
     770
    771771    int getlowbound(int iBoundNum = 0) const
    772772    {
    … …  
    774774    };
    775775
    776    
     776
    777777    int gethighbound(int iBoundNum = 0) const
    778778    {
    … …  
    788788    };
    789789
    790    
     790
    791791    const_raw_vector<T> getvector(int iStart, int iEnd) const
    792792    {
    … …  
    10341034namespace amp
    10351035{
    1036     class exception {}; 
     1036    class exception {};
    10371037    class incorrectPrecision    : public exception {};
    10381038    class overflow              : public exception {};
    … …  
    10431043    class internalError         : public exception {};
    10441044    class domainError           : public exception {};
    1045    
     1045
    10461046    typedef unsigned long unsigned32;
    10471047    typedef signed long   signed32;
    1048    
     1048
    10491049    struct mpfr_record
    10501050    {
    … …  
    10541054        mpfr_record *next;
    10551055    };
    1056    
     1056
    10571057    typedef mpfr_record* mpfr_record_ptr;
    1058    
     1058
    10591059    //
    10601060    // storage for mpfr_t instances
    … …  
    10621062    class mpfr_storage
    10631063    {
    1064     public:       
     1064    public:
    10651065        static mpfr_record* newMpfr(unsigned int Precision);
    10661066        static void deleteMpfr(mpfr_record* ref);
    … …  
    10701070        static mpfr_record_ptr& getList(unsigned int Precision);
    10711071    };
    1072    
     1072
    10731073    //
    10741074    // mpfr_t reference
    … …  
    10811081        mpfr_reference& operator= (const mpfr_reference &r);
    10821082        ~mpfr_reference();
    1083        
     1083
    10841084        void initialize(int Precision);
    10851085        void free();
    1086        
     1086
    10871087        mpfr_srcptr getReadPtr() const;
    10881088        mpfr_ptr getWritePtr();
    … …  
    10901090        mpfr_record *ref;
    10911091    };
    1092        
     1092
    10931093    //
    10941094    // ampf template
    … …  
    11071107                mpfr_storage::deleteMpfr(rval);
    11081108        }
    1109        
     1109
    11101110        //
    11111111        // Initializing
    … …  
    11131113        ampf ()                 { InitializeAsZero(); }
    11141114        ampf(mpfr_record *v)    { rval = v; }
    1115        
     1115
    11161116        ampf (long double v)    { InitializeAsDouble(v); }
    11171117        ampf (double v)         { InitializeAsDouble(v); }
    … …  
    11251125        ampf (signed char v)    { InitializeAsSLong(v); }
    11261126        ampf (unsigned char v)  { InitializeAsULong(v); }
    1127        
     1127
    11281128        //
    11291129        // initializing from string
    … …  
    11321132        ampf (const std::string &s) { InitializeAsString(s.c_str()); }
    11331133        ampf (const char *s)        { InitializeAsString(s); }
    1134        
     1134
    11351135        //
    11361136        // copy constructors
    … …  
    11921192        }
    11931193#endif
    1194        
     1194
    11951195        //
    11961196        // in-place operators
    … …  
    12011201        template<class T> ampf& operator*=(const T& v){ *this = *this * v; return *this; };
    12021202        template<class T> ampf& operator/=(const T& v){ *this = *this / v; return *this; };
    1203        
     1203
    12041204        //
    12051205        // MPFR access
    … …  
    12071207        mpfr_srcptr getReadPtr() const;
    12081208        mpfr_ptr getWritePtr();
    1209        
     1209
    12101210        //
    12111211        // properties and information
    … …  
    12241224        std::string toDec() const;
    12251225        char * toString() const;
    1226        
    1227        
     1226
     1227
    12281228        //
    12291229        // static methods
    1230         //       
     1230        //
    12311231        static const ampf getUlpOf(const ampf &x);
    12321232        static const ampf getUlp();
    … …  
    12461246        void InitializeAsDouble(long double v);
    12471247        void InitializeAsString(const char *s);
    1248        
     1248
    12491249        //mpfr_reference  ref;
    12501250        mpfr_record *rval;
    … …  
    12621262        if( Precision<32 )
    12631263            //throw incorrectPrecision();
    1264             WerrorS("incorrectPrecision");
     1264            WerrorS("incorrectPrecision");
    12651265    }
    12661266
    … …  
    12801280        mpfr_set_si(getWritePtr(), sv, GMP_RNDN);
    12811281    }
    1282    
     1282
    12831283    template<unsigned int Precision>
    12841284    void ampf<Precision>::InitializeAsULong(unsigned long v)
    … …  
    12881288        mpfr_set_ui(getWritePtr(), v, GMP_RNDN);
    12891289    }
    1290                            
     1290
    12911291    template<unsigned int Precision>
    12921292    void ampf<Precision>::InitializeAsDouble(long double v)
    … …  
    13031303        rval = mpfr_storage::newMpfr(Precision);
    13041304        mpfr_strtofr(getWritePtr(), s, NULL, 0, GMP_RNDN);
    1305     } 
     1305    }
    13061306
    13071307    template<unsigned int Precision>
    … …  
    13101310        return rval->value;
    13111311    }
    1312    
     1312
    13131313    template<unsigned int Precision>
    13141314    mpfr_ptr ampf<Precision>::getWritePtr()
    … …  
    13221322        return rval->value;
    13231323    }
    1324        
     1324
    13251325    template<unsigned int Precision>
    13261326    bool ampf<Precision>::isFiniteNumber() const
    … …  
    13361336        return mpfr_sgn(getReadPtr())>0;
    13371337    }
    1338    
     1338
    13391339    template<unsigned int Precision>
    13401340    bool ampf<Precision>::isZero() const
    … …  
    13421342        return mpfr_zero_p(getReadPtr())!=0;
    13431343    }
    1344    
     1344
    13451345    template<unsigned int Precision>
    13461346    bool ampf<Precision>::isNegativeNumber() const
    … …  
    13561356        return getUlpOf(*this);
    13571357    }
    1358    
     1358
    13591359    template<unsigned int Precision>
    13601360    double ampf<Precision>::toDouble() const
    … …  
    13621362        return mpfr_get_d(getReadPtr(), GMP_RNDN);
    13631363    }
    1364    
     1364
    13651365    template<unsigned int Precision>
    13661366    std::string ampf<Precision>::toHex() const
    … …  
    13791379            return r;
    13801380        }
    1381        
     1381
    13821382        //
    13831383        // general case
    … …  
    13931393        iexpval = expval;
    13941394        if( iexpval!=expval )
    1395             throw internalError();
     1395        //    throw internalError();
     1396            WerrorS("internalError");
    13961397        sprintf(buf_e, "%ld", long(iexpval));
    13971398        if( *ptr=='-' )
    … …  
    14071408        return r;
    14081409    }
    1409    
     1410
    14101411    template<unsigned int Precision>
    14111412    std::string ampf<Precision>::toDec() const
    14121413    {
    14131414        // TODO: advanced output formatting (zero, integers)
    1414        
     1415
    14151416        //
    14161417        // some special cases
    … …  
    14261427            return r;
    14271428        }
    1428        
     1429
    14291430        //
    14301431        // general case
    … …  
    14401441        iexpval = expval;
    14411442        if( iexpval!=expval )
    1442             throw internalError();
     1443        //    throw internalError();
     1444            WerrorS("internalError");
    14431445        sprintf(buf_e, "%ld", long(iexpval));
    14441446        if( *ptr=='-' )
    … …  
    14701472            return toString_Block;
    14711473        }
    1472        
     1474
    14731475        //
    14741476        // general case
    … …  
    14851487        if( iexpval!=expval )
    14861488            //throw internalError();
    1487             WerrorS("internalError");
     1489            WerrorS("internalError");
    14881490        sprintf(buf_e, "%ld", long(iexpval));
    14891491        if( *ptr=='-' )
    … …  
    15201522        return r;
    15211523    }
    1522        
     1524
    15231525    template<unsigned int Precision>
    15241526    const ampf<Precision> ampf<Precision>::getUlp()
    … …  
    15291531        return r;
    15301532    }
    1531    
     1533
    15321534    template<unsigned int Precision>
    15331535    const ampf<Precision> ampf<Precision>::getUlp256()
    … …  
    15431545        return r;
    15441546    }
    1545    
     1547
    15461548    template<unsigned int Precision>
    15471549    const ampf<Precision> ampf<Precision>::getUlp512()
    … …  
    15561558            GMP_RNDN);
    15571559        return r;
    1558     } 
     1560    }
    15591561
    15601562    template<unsigned int Precision>
    … …  
    15661568        return r;
    15671569    }
    1568    
     1570
    15691571    template<unsigned int Precision>
    15701572    const ampf<Precision> ampf<Precision>::getMinNumber()
    … …  
    15801582        return getUlp256();
    15811583    }
    1582    
     1584
    15831585    template<unsigned int Precision>
    15841586    const ampf<Precision> ampf<Precision>::getAlgoPascalMaxNumber()
    … …  
    15911593        return r;
    15921594    }
    1593    
     1595
    15941596    template<unsigned int Precision>
    15951597    const ampf<Precision> ampf<Precision>::getAlgoPascalMinNumber()
    … …  
    16101612        return r;
    16111613    }
    1612    
     1614
    16131615    //
    16141616    // comparison operators
    … …  
    16251627        return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())!=0;
    16261628    }
    1627    
     1629
    16281630    template<unsigned int Precision>
    16291631    const bool operator<(const ampf<Precision>& op1, const ampf<Precision>& op2)
    … …  
    16311633        return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())<0;
    16321634    }
    1633    
     1635
    16341636    template<unsigned int Precision>
    16351637    const bool operator>(const ampf<Precision>& op1, const ampf<Precision>& op2)
    … …  
    16371639        return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())>0;
    16381640    }
    1639    
     1641
    16401642    template<unsigned int Precision>
    16411643    const bool operator<=(const ampf<Precision>& op1, const ampf<Precision>& op2)
    … …  
    16431645        return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())<=0;
    16441646    }
    1645    
     1647
    16461648    template<unsigned int Precision>
    16471649    const bool operator>=(const ampf<Precision>& op1, const ampf<Precision>& op2)
    … …  
    16491651        return mpfr_cmp(op1.getReadPtr(), op2.getReadPtr())>=0;
    16501652    }
    1651    
     1653
    16521654    //
    16531655    // arithmetic operators
    … …  
    16581660        return op1;
    16591661    }
    1660    
     1662
    16611663    template<unsigned int Precision>
    16621664    const ampf<Precision> operator-(const ampf<Precision>& op1)
    … …  
    16661668        return v;
    16671669    }
    1668    
     1670
    16691671    template<unsigned int Precision>
    16701672    const ampf<Precision> operator+(const ampf<Precision>& op1, const ampf<Precision>& op2)
    … …  
    16741676        return v;
    16751677    }
    1676    
     1678
    16771679    template<unsigned int Precision>
    16781680    const ampf<Precision> operator-(const ampf<Precision>& op1, const ampf<Precision>& op2)
    … …  
    16821684        return v;
    16831685    }
    1684        
    1685    
     1686
     1687
    16861688    template<unsigned int Precision>
    16871689    const ampf<Precision> operator*(const ampf<Precision>& op1, const ampf<Precision>& op2)
    … …  
    16911693        return v;
    16921694    }
    1693    
     1695
    16941696    template<unsigned int Precision>
    16951697    const ampf<Precision> operator/(const ampf<Precision>& op1, const ampf<Precision>& op2)
    … …  
    16991701        return v;
    17001702    }
    1701  
     1703
    17021704    //
    17031705    // basic functions
    … …  
    17661768        mpfr_trunc(tmp.getWritePtr(), x.getReadPtr());
    17671769        if( mpfr_integer_p(tmp.getReadPtr())==0 )
    1768             throw invalidConversion();
     1770            //throw invalidConversion();
     1771            WerrorS("internalError");
    17691772        mpfr_clear_erangeflag();
    17701773        r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
    17711774        if( mpfr_erangeflag_p()!=0 )
    1772             throw invalidConversion();
     1775            //throw invalidConversion();
     1776            WerrorS("internalError");
    17731777        return r;
    17741778    }
    … …  
    17901794        mpfr_floor(tmp.getWritePtr(), x.getReadPtr());
    17911795        if( mpfr_integer_p(tmp.getReadPtr())==0 )
    1792             throw invalidConversion();
     1796            //throw invalidConversion();
     1797            WerrorS("internalError");
    17931798        mpfr_clear_erangeflag();
    17941799        r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
    17951800        if( mpfr_erangeflag_p()!=0 )
    1796             throw invalidConversion();
     1801            //throw invalidConversion();
     1802            WerrorS("internalError");
    17971803        return r;
    17981804    }
    … …  
    18051811        mpfr_ceil(tmp.getWritePtr(), x.getReadPtr());
    18061812        if( mpfr_integer_p(tmp.getReadPtr())==0 )
    1807             throw invalidConversion();
     1813            //throw invalidConversion();
     1814            WerrorS("internalError");
    18081815        mpfr_clear_erangeflag();
    18091816        r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
    18101817        if( mpfr_erangeflag_p()!=0 )
    1811             throw invalidConversion();
     1818            //throw invalidConversion();
     1819            WerrorS("internalError");
    18121820        return r;
    18131821    }
    … …  
    18201828        mpfr_round(tmp.getWritePtr(), x.getReadPtr());
    18211829        if( mpfr_integer_p(tmp.getReadPtr())==0 )
    1822             throw invalidConversion();
     1830            //throw invalidConversion();
     1831            WerrorS("internalError");
    18231832        mpfr_clear_erangeflag();
    18241833        r = mpfr_get_si(tmp.getReadPtr(), GMP_RNDN);
    18251834        if( mpfr_erangeflag_p()!=0 )
    1826             throw invalidConversion();
     1835            //throw invalidConversion();
     1836            WerrorS("internalError");
    18271837        return r;
    18281838    }
    … …  
    18341844        ampf<Precision> r;
    18351845        if( !x.isFiniteNumber() )
    1836             throw invalidConversion();
     1846            //throw invalidConversion();
     1847            WerrorS("internalError");
    18371848        if( x.isZero() )
    18381849        {
    … …  
    18551866        return r;
    18561867    }
    1857  
     1868
    18581869    //
    18591870    // different types of arguments
    … …  
    19301941    __AMP_BINARY_OPF(long double)
    19311942    #undef __AMP_BINARY_OPF
    1932    
     1943
    19331944    //
    19341945    // transcendent functions
    … …  
    20152026        return v;
    20162027    }
    2017    
     2028
    20182029    template<unsigned int Precision>
    20192030    const ampf<Precision> log(const ampf<Precision> &x)
    … …  
    20892100        campf():x(0),y(0){};
    20902101        campf(long double v)    { x=v; y=0; }
    2091         campf(double v)         { x=v; y=0; } 
     2102        campf(double v)         { x=v; y=0; }
    20922103        campf(float v)          { x=v; y=0; }
    20932104        campf(signed long v)    { x=v; y=0; }
    … …  
    21382149        ampf<Precision> x, y;
    21392150    };
    2140    
     2151
    21412152    //
    21422153    // complex operations
    … …  
    21452156    const bool operator==(const campf<Precision>& lhs, const campf<Precision>& rhs)
    21462157    { return lhs.x==rhs.x && lhs.y==rhs.y; }
    2147    
     2158
    21482159    template<unsigned int Precision>
    21492160    const bool operator!=(const campf<Precision>& lhs, const campf<Precision>& rhs)
    21502161    { return lhs.x!=rhs.x || lhs.y!=rhs.y; }
    2151    
     2162
    21522163    template<unsigned int Precision>
    21532164    const campf<Precision> operator+(const campf<Precision>& lhs)
    21542165    { return lhs; }
    2155    
     2166
    21562167    template<unsigned int Precision>
    21572168    campf<Precision>& operator+=(campf<Precision>& lhs, const campf<Precision>& rhs)
    21582169    { lhs.x += rhs.x; lhs.y += rhs.y; return lhs; }
    2159    
     2170
    21602171    template<unsigned int Precision>
    21612172    const campf<Precision> operator+(const campf<Precision>& lhs, const campf<Precision>& rhs)
    21622173    { campf<Precision> r = lhs; r += rhs; return r; }
    2163    
     2174
    21642175    template<unsigned int Precision>
    21652176    const campf<Precision> operator-(const campf<Precision>& lhs)
    21662177    { return campf<Precision>(-lhs.x, -lhs.y); }
    2167    
     2178
    21682179    template<unsigned int Precision>
    21692180    campf<Precision>& operator-=(campf<Precision>& lhs, const campf<Precision>& rhs)
    21702181    { lhs.x -= rhs.x; lhs.y -= rhs.y; return lhs; }
    2171    
     2182
    21722183    template<unsigned int Precision>
    21732184    const campf<Precision> operator-(const campf<Precision>& lhs, const campf<Precision>& rhs)
    21742185    { campf<Precision> r = lhs; r -= rhs; return r; }
    2175    
     2186
    21762187    template<unsigned int Precision>
    21772188    campf<Precision>& operator*=(campf<Precision>& lhs, const campf<Precision>& rhs)
    … …  
    21822193        return lhs;
    21832194    }
    2184    
     2195
    21852196    template<unsigned int Precision>
    21862197    const campf<Precision> operator*(const campf<Precision>& lhs, const campf<Precision>& rhs)
    21872198    { campf<Precision> r = lhs; r *= rhs; return r; }
    2188    
     2199
    21892200    template<unsigned int Precision>
    21902201    const campf<Precision> operator/(const campf<Precision>& lhs, const campf<Precision>& rhs)
    … …  
    22092220        return result;
    22102221    }
    2211    
     2222
    22122223    template<unsigned int Precision>
    22132224    campf<Precision>& operator/=(campf<Precision>& lhs, const campf<Precision>& rhs)
    … …  
    22162227        return lhs;
    22172228    }
    2218    
     2229
    22192230    template<unsigned int Precision>
    22202231    const ampf<Precision> abscomplex(const campf<Precision> &z)
    22212232    {
    22222233        ampf<Precision> w, xabs, yabs, v;
    2223    
     2234
    22242235        xabs = abs(z.x);
    22252236        yabs = abs(z.y);
    22262237        w = xabs>yabs ? xabs : yabs;
    2227         v = xabs<yabs ? xabs : yabs; 
     2238        v = xabs<yabs ? xabs : yabs;
    22282239        if( v==0 )
    22292240            return w;
    … …  
    22342245        }
    22352246    }
    2236    
     2247
    22372248    template<unsigned int Precision>
    22382249    const campf<Precision> conj(const campf<Precision> &z)
    22392250    {
    2240         return campf<Precision>(z.x, -z.y); 
    2241     }
    2242    
     2251        return campf<Precision>(z.x, -z.y);
     2252    }
     2253
    22432254    template<unsigned int Precision>
    22442255    const campf<Precision> csqr(const campf<Precision> &z)
    22452256    {
    2246         ampf<Precision> t = z.x*z.y; 
    2247         return campf<Precision>(sqr(z.x)-sqr(z.y), t+t); 
    2248     }
    2249    
     2257        ampf<Precision> t = z.x*z.y;
     2258        return campf<Precision>(sqr(z.x)-sqr(z.y), t+t);
     2259    }
     2260
    22502261    //
    22512262    // different types of arguments
    … …  
    22932304        template<unsigned int Precision>                   bool operator==(const campf<Precision>& op1, const type& op2)             { return op1.x==op2 && op1.y==0; }   \
    22942305        template<unsigned int Precision>                   bool operator!=(const type& op1,             const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; }   \
    2295         template<unsigned int Precision>                   bool operator!=(const campf<Precision>& op1, const type& op2)             { return op1.x!=op2 || op1.y!=0; }   
     2306        template<unsigned int Precision>                   bool operator!=(const campf<Precision>& op1, const type& op2)             { return op1.x!=op2 || op1.y!=0; }
    22962307    __AMP_BINARY_OPF(float)
    22972308    __AMP_BINARY_OPF(double)
    … …  
    22992310    __AMP_BINARY_OPF(ampf<Precision>)
    23002311    #undef __AMP_BINARY_OPF
    2301    
     2312
    23022313    //
    23032314    // Real linear algebra
    … …  
    23272338            return r;
    23282339        }
    2329         catch(...)
    2330         {
    2331             if( r!=NULL )
    2332                 mpfr_storage::deleteMpfr(r);
    2333             if( t!=NULL )
    2334                 mpfr_storage::deleteMpfr(t);
    2335             throw;
    2336         }
     2340        //catch(...)
     2341        //{
     2342        //    if( r!=NULL )
     2343        //        mpfr_storage::deleteMpfr(r);
     2344        //    if( t!=NULL )
     2345        //        mpfr_storage::deleteMpfr(t);
     2346        //    throw;
     2347        //}
    23372348    }
    23382349
    … …  
    23532364        }
    23542365    }
    2355    
     2366
    23562367    template<unsigned int Precision>
    23572368    void vMoveNeg(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
    … …  
    23702381        }
    23712382    }
    2372    
     2383
    23732384    template<unsigned int Precision, class T2>
    23742385    void vMove(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc, T2 alpha)
    … …  
    23882399        }
    23892400    }
    2390    
     2401
    23912402    template<unsigned int Precision>
    23922403    void vAdd(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
    … …  
    24052416        }
    24062417    }
    2407    
     2418
    24082419    template<unsigned int Precision, class T2>
    24092420    void vAdd(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc, T2 alpha)
    … …  
    24242435        }
    24252436    }
    2426    
     2437
    24272438    template<unsigned int Precision>
    24282439    void vSub(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc)
    … …  
    24412452        }
    24422453    }
    2443    
     2454
    24442455    template<unsigned int Precision, class T2>
    24452456    void vSub(ap::raw_vector< ampf<Precision> > vDst, ap::const_raw_vector< ampf<Precision> > vSrc, T2 alpha)
    … …  
    24472458        vAdd(vDst, vSrc, -alpha);
    24482459    }
    2449    
     2460
    24502461    template<unsigned int Precision, class T2>
    24512462    void vMul(ap::raw_vector< ampf<Precision> > vDst, T2 alpha)
    … …  
    24622473    }
    24632474}
    2464  
     2475
    24652476#endif
    24662477
    … …  
    25872598
    25882599
    2589        
     2600
    25902601        //
    25912602        // Executable Statements ..
    … …  
    25962607            return;
    25972608        }
    2598        
     2609
    25992610        //
    26002611        // XNORM = DNRM2( N-1, X, INCX )
    … …  
    26172628        if( xnorm==0 )
    26182629        {
    2619            
     2630
    26202631            //
    26212632            // H  =  I
    … …  
    26242635            return;
    26252636        }
    2626        
     2637
    26272638        //
    26282639        // general case
    … …  
    26882699            return;
    26892700        }
    2690        
     2701
    26912702        //
    26922703        // w := C' * v
    … …  
    27022713            ap::vadd(work.getvector(n1, n2), c.getrow(i, n1, n2), t);
    27032714        }
    2704        
     2715
    27052716        //
    27062717        // C := C - tau * v * w'
    … …  
    27612772            return;
    27622773        }
    2763        
     2774
    27642775        //
    27652776        // w := C * v
    … …  
    27712782            work(i) = t;
    27722783        }
    2773        
     2784
    27742785        //
    27752786        // C := C - w * v'
    … …  
    29923003
    29933004
    2994        
     3005
    29953006        //
    29963007        // Prepare
    … …  
    30163027        if( m>=n )
    30173028        {
    3018            
     3029
    30193030            //
    30203031            // Reduce to upper bidiagonal form
    … …  
    30223033            for(i=0; i<=n-1; i++)
    30233034            {
    3024                
     3035
    30253036                //
    30263037                // Generate elementary reflector H(i) to annihilate A(i+1:m-1,i)
    … …  
    30313042                ap::vmove(a.getcolumn(i, i, m-1), t.getvector(1, m-i));
    30323043                t(1) = 1;
    3033                
     3044
    30343045                //
    30353046                // Apply H(i) to A(i:m-1,i+1:n-1) from the left
    … …  
    30383049                if( i<n-1 )
    30393050                {
    3040                    
     3051
    30413052                    //
    30423053                    // Generate elementary reflector G(i) to annihilate
    … …  
    30483059                    ap::vmove(a.getrow(i, i+1, n-1), t.getvector(1, n-1-i));
    30493060                    t(1) = 1;
    3050                    
     3061
    30513062                    //
    30523063                    // Apply G(i) to A(i+1:m-1,i+1:n-1) from the right
    … …  
    30623073        else
    30633074        {
    3064            
     3075
    30653076            //
    30663077            // Reduce to lower bidiagonal form
    … …  
    30683079            for(i=0; i<=m-1; i++)
    30693080            {
    3070                
     3081
    30713082                //
    30723083                // Generate elementary reflector G(i) to annihilate A(i,i+1:n-1)
    … …  
    30773088                ap::vmove(a.getrow(i, i, n-1), t.getvector(1, n-i));
    30783089                t(1) = 1;
    3079                
     3090
    30803091                //
    30813092                // Apply G(i) to A(i+1:m-1,i:n-1) from the right
    … …  
    30843095                if( i<m-1 )
    30853096                {
    3086                    
     3097
    30873098                    //
    30883099                    // Generate elementary reflector H(i) to annihilate
    … …  
    30943105                    ap::vmove(a.getcolumn(i, i+1, m-1), t.getvector(1, m-1-i));
    30953106                    t(1) = 1;
    3096                    
     3107
    30973108                    //
    30983109                    // Apply H(i) to A(i+1:m-1,i+1:n-1) from the left
    … …  
    31483159            return;
    31493160        }
    3150        
     3161
    31513162        //
    31523163        // prepare Q
    … …  
    31673178            }
    31683179        }
    3169        
     3180
    31703181        //
    31713182        // Calculate
    … …  
    32293240        }
    32303241        ap::ap_error::make_assertion(fromtheright && zcolumns==m || !fromtheright && zrows==m);
    3231        
     3242
    32323243        //
    32333244        // init
    … …  
    32403251        if( m>=n )
    32413252        {
    3242            
     3253
    32433254            //
    32443255            // setup
    … …  
    32633274                istep = -istep;
    32643275            }
    3265            
     3276
    32663277            //
    32673278            // Process
    … …  
    32863297        else
    32873298        {
    3288            
     3299
    32893300            //
    32903301            // setup
    … …  
    33093320                istep = -istep;
    33103321            }
    3311            
     3322
    33123323            //
    33133324            // Process
    … …  
    33753386            return;
    33763387        }
    3377        
     3388
    33783389        //
    33793390        // prepare PT
    … …  
    33943405            }
    33953406        }
    3396        
     3407
    33973408        //
    33983409        // Calculate
    … …  
    34563467        }
    34573468        ap::ap_error::make_assertion(fromtheright && zcolumns==n || !fromtheright && zrows==n);
    3458        
     3469
    34593470        //
    34603471        // init
    … …  
    34693480        if( m>=n )
    34703481        {
    3471            
     3482
    34723483            //
    34733484            // setup
    … …  
    34923503                istep = -istep;
    34933504            }
    3494            
     3505
    34953506            //
    34963507            // Process
    … …  
    35183529        else
    35193530        {
    3520            
     3531
    35213532            //
    35223533            // setup
    … …  
    35413552                istep = -istep;
    35423553            }
    3543            
     3554
    35443555            //
    35453556            // Process
    … …  
    36603671        if( m>=n )
    36613672        {
    3662            
     3673
    36633674            //
    36643675            // Reduce to upper bidiagonal form
    … …  
    36663677            for(i=1; i<=n; i++)
    36673678            {
    3668                
     3679
    36693680                //
    36703681                // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
    … …  
    36763687                ap::vmove(a.getcolumn(i, i, m), t.getvector(1, mmip1));
    36773688                t(1) = 1;
    3678                
     3689
    36793690                //
    36803691                // Apply H(i) to A(i:m,i+1:n) from the left
    … …  
    36833694                if( i<n )
    36843695                {
    3685                    
     3696
    36863697                    //
    36873698                    // Generate elementary reflector G(i) to annihilate
    … …  
    36953706                    ap::vmove(a.getrow(i, ip1, n), t.getvector(1, nmi));
    36963707                    t(1) = 1;
    3697                    
     3708
    36983709                    //
    36993710                    // Apply G(i) to A(i+1:m,i+1:n) from the right
    … …  
    37093720        else
    37103721        {
    3711            
     3722
    37123723            //
    37133724            // Reduce to lower bidiagonal form
    … …  
    37153726            for(i=1; i<=m; i++)
    37163727            {
    3717                
     3728
    37183729                //
    37193730                // Generate elementary reflector G(i) to annihilate A(i,i+1:n)
    … …  
    37253736                ap::vmove(a.getrow(i, i, n), t.getvector(1, nmip1));
    37263737                t(1) = 1;
    3727                
     3738
    37283739                //
    37293740                // Apply G(i) to A(i+1:m,i:n) from the right
    … …  
    37323743                if( i<m )
    37333744                {
    3734                    
     3745
    37353746                    //
    37363747                    // Generate elementary reflector H(i) to annihilate
    … …  
    37443755                    ap::vmove(a.getcolumn(i, ip1, m), t.getvector(1, mmi));
    37453756                    t(1) = 1;
    3746                    
     3757
    37473758                    //
    37483759                    // Apply H(i) to A(i+1:m,i+1:n) from the left
    … …  
    37843795            return;
    37853796        }
    3786        
     3797
    37873798        //
    37883799        // init
    … …  
    37913802        v.setbounds(1, m);
    37923803        work.setbounds(1, qcolumns);
    3793        
     3804
    37943805        //
    37953806        // prepare Q
    … …  
    38643875        }
    38653876        ap::ap_error::make_assertion(fromtheright && zcolumns==m || !fromtheright && zrows==m);
    3866        
     3877
    38673878        //
    38683879        // init
    … …  
    38753886        if( m>=n )
    38763887        {
    3877            
     3888
    38783889            //
    38793890            // setup
    … …  
    38983909                istep = -istep;
    38993910            }
    3900            
     3911
    39013912            //
    39023913            // Process
    … …  
    39223933        else
    39233934        {
    3924            
     3935
    39253936            //
    39263937            // setup
    … …  
    39453956                istep = -istep;
    39463957            }
    3947            
     3958
    39483959            //
    39493960            // Process
    … …  
    39994010            return;
    40004011        }
    4001        
     4012
    40024013        //
    40034014        // init
    … …  
    40064017        v.setbounds(1, n);
    40074018        work.setbounds(1, ptrows);
    4008        
     4019
    40094020        //
    40104021        // prepare PT
    … …  
    40794090        }
    40804091        ap::ap_error::make_assertion(fromtheright && zcolumns==n || !fromtheright && zrows==n);
    4081        
     4092
    40824093        //
    40834094        // init
    … …  
    40924103        if( m>=n )
    40934104        {
    4094            
     4105
    40954106            //
    40964107            // setup
    … …  
    41154126                istep = -istep;
    41164127            }
    4117            
     4128
    41184129            //
    41194130            // Process
    … …  
    41434154        else
    41444155        {
    4145            
     4156
    41464157            //
    41474158            // setup
    … …  
    41664177                istep = -istep;
    41674178            }
    4168            
     4179
    41694180            //
    41704181            // Process
    … …  
    43804391        t.setbounds(1, m);
    43814392        tau.setbounds(0, minmn-1);
    4382        
     4393
    43834394        //
    43844395        // Test the input arguments
    … …  
    43874398        for(i=0; i<=k-1; i++)
    43884399        {
    4389            
     4400
    43904401            //
    43914402            // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
    … …  
    43984409            if( i<n )
    43994410            {
    4400                
     4411
    44014412                //
    44024413                // Apply H(i) to A(i:m-1,i+1:n-1) from the left
    … …  
    44494460            return;
    44504461        }
    4451        
     4462
    44524463        //
    44534464        // init
    … …  
    44724483            }
    44734484        }
    4474        
     4485
    44754486        //
    44764487        // unpack Q
    … …  
    44784489        for(i=k-1; i>=0; i--)
    44794490        {
    4480            
     4491
    44814492            //
    44824493            // Apply H(i)
    … …  
    45574568        t.setbounds(1, m);
    45584569        tau.setbounds(1, minmn);
    4559        
     4570
    45604571        //
    45614572        // Test the input arguments
    … …  
    45644575        for(i=1; i<=k; i++)
    45654576        {
    4566            
     4577
    45674578            //
    45684579            // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
    … …  
    45764587            if( i<n )
    45774588            {
    4578                
     4589
    45794590                //
    45804591                // Apply H(i) to A(i:m,i+1:n) from the left
    … …  
    46114622            return;
    46124623        }
    4613        
     4624
    46144625        //
    46154626        // init
    … …  
    46344645            }
    46354646        }
    4636        
     4647
    46374648        //
    46384649        // unpack Q
    … …  
    46404651        for(i=k; i>=1; i--)
    46414652        {
    4642            
     4653
    46434654            //
    46444655            // Apply H(i)
    … …  
    46784689        q.setbounds(1, m, 1, m);
    46794690        r.setbounds(1, m, 1, n);
    4680        
     4691
    46814692        //
    46824693        // QRDecomposition
    46834694        //
    46844695        qrdecomposition<Precision>(a, m, n, tau);
    4685        
     4696
    46864697        //
    46874698        // R
    … …  
    46994710            ap::vmove(r.getrow(i, i, n), a.getrow(i, i, n));
    47004711        }
    4701        
     4712
    47024713        //
    47034714        // Q
    … …  
    48424853        for(i=0; i<=k-1; i++)
    48434854        {
    4844            
     4855
    48454856            //
    48464857            // Generate elementary reflector H(i) to annihilate A(i,i+1:n-1)
    … …  
    48534864            if( i<n )
    48544865            {
    4855                
     4866
    48564867                //
    48574868                // Apply H(i) to A(i+1:m,i:n) from the right
    … …  
    49044915            return;
    49054916        }
    4906        
     4917
    49074918        //
    49084919        // init
    … …  
    49274938            }
    49284939        }
    4929        
     4940
    49304941        //
    49314942        // unpack Q
    … …  
    49334944        for(i=k-1; i>=0; i--)
    49344945        {
    4935            
     4946
    49364947            //
    49374948            // Apply H(i)
    … …  
    50155026        t.setbounds(1, n);
    50165027        tau.setbounds(1, minmn);
    5017        
     5028
    50185029        //
    50195030        // Test the input arguments
    … …  
    50225033        for(i=1; i<=k; i++)
    50235034        {
    5024            
     5035
    50255036            //
    50265037            // Generate elementary reflector H(i) to annihilate A(i,i+1:n)
    … …  
    50345045            if( i<n )
    50355046            {
    5036                
     5047
    50375048                //
    50385049                // Apply H(i) to A(i+1:m,i:n) from the right
    … …  
    50705081            return;
    50715082        }
    5072        
     5083
    50735084        //
    50745085        // init
    … …  
    50935104            }
    50945105        }
    5095        
     5106
    50965107        //
    50975108        // unpack Q
    … …  
    50995110        for(i=k; i>=1; i--)
    51005111        {
    5101            
     5112
    51025113            //
    51035114            // Apply H(i)
    … …  
    51325143        q.setbounds(1, n, 1, n);
    51335144        l.setbounds(1, m, 1, n);
    5134        
     5145
    51355146        //
    51365147        // LQDecomposition
    51375148        //
    51385149        lqdecomposition<Precision>(a, m, n, tau);
    5139        
     5150
    51405151        //
    51415152        // L
    … …  
    51555166            }
    51565167        }
    5157        
     5168
    51585169        //
    51595170        // Q
    … …  
    55665577        if( !trans )
    55675578        {
    5568            
     5579
    55695580            //
    55705581            // y := alpha*A*x + beta*y;
    … …  
    55765587            ap::ap_error::make_assertion(j2-j1==ix2-ix1);
    55775588            ap::ap_error::make_assertion(i2-i1==iy2-iy1);
    5578            
     5589
    55795590            //
    55805591            // beta*y
    … …  
    55915602                ap::vmul(y.getvector(iy1, iy2), beta);
    55925603            }
    5593            
     5604
    55945605            //
    55955606            // alpha*A*x
    … …  
    56035614        else
    56045615        {
    5605            
     5616
    56065617            //
    56075618            // y := alpha*A'*x + beta*y;
    … …  
    56135624            ap::ap_error::make_assertion(i2-i1==ix2-ix1);
    56145625            ap::ap_error::make_assertion(j2-j1==iy2-iy1);
    5615            
     5626
    56165627            //
    56175628            // beta*y
    … …  
    56285639                ap::vmul(y.getvector(iy1, iy2), beta);
    56295640            }
    5630            
     5641
    56315642            //
    56325643            // alpha*A'*x
    … …  
    57045715
    57055716
    5706        
     5717
    57075718        //
    57085719        // Setup
    … …  
    57355746        crows = arows;
    57365747        ccols = bcols;
    5737        
     5748
    57385749        //
    57395750        // Test WORK
    … …  
    57445755        work(1) = 0;
    57455756        work(i) = 0;
    5746        
     5757
    57475758        //
    57485759        // Prepare C
    … …  
    57655776            }
    57665777        }
    5767        
     5778
    57685779        //
    57695780        // A*B
    … …  
    57825793            return;
    57835794        }
    5784        
     5795
    57855796        //
    57865797        // A*B'
    … …  
    58135824            }
    58145825        }
    5815        
     5826
    58165827        //
    58175828        // A'*B
    … …  
    58305841            return;
    58315842        }
    5832        
     5843
    58335844        //
    58345845        // A'*B'
    … …  
    59946005            return;
    59956006        }
    5996        
     6007
    59976008        //
    59986009        // Form  P * A
    … …  
    60026013            if( n1!=n2 )
    60036014            {
    6004                
     6015
    60056016                //
    60066017                // Common case: N1<>N2
    … …  
    60236034            else
    60246035            {
    6025                
     6036
    60266037                //
    60276038                // Special case: N1=N2
    … …  
    60446055            if( n1!=n2 )
    60456056            {
    6046                
     6057
    60476058                //
    60486059                // Common case: N1<>N2
    … …  
    60656076            else
    60666077            {
    6067                
     6078
    60686079                //
    60696080                // Special case: N1=N2
    … …  
    61286139
    61296140
    6130        
     6141
    61316142        //
    61326143        // Form A * P'
    … …  
    61366147            if( m1!=m2 )
    61376148            {
    6138                
     6149
    61396150                //
    61406151                // Common case: M1<>M2
    … …  
    61576168            else
    61586169            {
    6159                
     6170
    61606171                //
    61616172                // Special case: M1=M2
    … …  
    61786189            if( m1!=m2 )
    61796190            {
    6180                
     6191
    61816192                //
    61826193                // Common case: M1<>M2
    … …  
    61996210            else
    62006211            {
    6201                
     6212
    62026213                //
    62036214                // Special case: M1=M2
    … …  
    66186629            return result;
    66196630        }
    6620        
     6631
    66216632        //
    66226633        // init
    … …  
    66356646        rightside = true;
    66366647        fwddir = true;
    6637        
     6648
    66386649        //
    66396650        // resize E from N-1 to N
    … …  
    66516662        e(n) = 0;
    66526663        idir = 0;
    6653        
     6664
    66546665        //
    66556666        // Get machine constants
    … …  
    66576668        eps = amp::ampf<Precision>::getAlgoPascalEpsilon();
    66586669        unfl = amp::ampf<Precision>::getAlgoPascalMinNumber();
    6659        
     6670
    66606671        //
    66616672        // If matrix lower bidiagonal, rotate to be upper bidiagonal
    … …  
    66736684                work1(i) = sn;
    66746685            }
    6675            
     6686
    66766687            //
    66776688            // Update singular vectors if desired
    … …  
    66866697            }
    66876698        }
    6688        
     6699
    66896700        //
    66906701        // Compute singular values to relative accuracy TOL
    … …  
    66986709            tol = -tol;
    66996710        }
    6700        
     6711
    67016712        //
    67026713        // Compute approximate maximum, minimum singular values
    … …  
    67146725        if( tol>=0 )
    67156726        {
    6716            
     6727
    67176728            //
    67186729            // Relative accuracy desired
    … …  
    67376748        else
    67386749        {
    6739            
     6750
    67406751            //
    67416752            // Absolute accuracy desired
    … …  
    67436754            thresh = amp::maximum<Precision>(amp::abs<Precision>(tol)*smax, maxitr*n*n*unfl);
    67446755        }
    6745        
     6756
    67466757        //
    67476758        // Prepare for main iteration loop for the singular values
    … …  
    67536764        oldll = -1;
    67546765        oldm = -1;
    6755        
     6766
    67566767        //
    67576768        // M points to last element of unconverged part of matrix
    67586769        //
    67596770        m = n;
    6760        
     6771
    67616772        //
    67626773        // Begin main iteration loop
    … …  
    67646775        while( true )
    67656776        {
    6766            
     6777
    67676778            //
    67686779            // Check for convergence or exceeding iteration count
    … …  
    67776788                return result;
    67786789            }
    6779            
     6790
    67806791            //
    67816792            // Find diagonal block of matrix to work on
    … …  
    68116822            else
    68126823            {
    6813                
     6824
    68146825                //
    68156826                // Matrix splits since E(LL) = 0
    … …  
    68186829                if( ll==m-1 )
    68196830                {
    6820                    
     6831
    68216832                    //
    68226833                    // Convergence of bottom singular value, return to top of loop
    … …  
    68276838            }
    68286839            ll = ll+1;
    6829            
     6840
    68306841            //
    68316842            // E(LL) through E(M-1) are nonzero, E(LL-1) is zero
    … …  
    68336844            if( ll==m-1 )
    68346845            {
    6835                
     6846
    68366847                //
    68376848                // 2 by 2 block, handle separately
    … …  
    68416852                e(m-1) = 0;
    68426853                d(m) = sigmn;
    6843                
     6854
    68446855                //
    68456856                // Compute singular vectors, if desired
    … …  
    68786889                continue;
    68796890            }
    6880            
     6891
    68816892            //
    68826893            // If working on new submatrix, choose shift direction
    … …  
    69016912                if( amp::abs<Precision>(d(ll))>=amp::abs<Precision>(d(m)) )
    69026913                {
    6903                    
     6914
    69046915                    //
    69056916                    // Chase bulge from top (big end) to bottom (small end)
    … …  
    69096920                else
    69106921                {
    6911                    
     6922
    69126923                    //
    69136924                    // Chase bulge from bottom (big end) to top (small end)
    … …  
    69166927                }
    69176928            }
    6918            
     6929
    69196930            //
    69206931            // Apply convergence tests
    … …  
    69226933            if( idir==1 )
    69236934            {
    6924                
     6935
    69256936                //
    69266937                // Run convergence test in forward direction
    … …  
    69346945                if( tol>=0 )
    69356946                {
    6936                    
     6947
    69376948                    //
    69386949                    // If relative accuracy desired,
    … …  
    69626973            else
    69636974            {
    6964                
     6975
    69656976                //
    69666977                // Run convergence test in backward direction
    … …  
    69746985                if( tol>=0 )
    69756986                {
    6976                    
     6987
    69776988                    //
    69786989                    // If relative accuracy desired,
    … …  
    70027013            oldll = ll;
    70037014            oldm = m;
    7004            
     7015
    70057016            //
    70067017            // Compute shift.  First, test if shifting would ruin relative
    … …  
    70097020            if( tol>=0 && n*tol*(sminl/smax)<=amp::maximum<Precision>(eps, amp::ampf<Precision>("0.01")*tol) )
    70107021            {
    7011                
     7022
    70127023                //
    70137024                // Use a zero shift to avoid loss of relative accuracy
    … …  
    70177028            else
    70187029            {
    7019                
     7030
    70207031                //
    70217032                // Compute the shift from 2-by-2 block at end of matrix
    … …  
    70317042                    svd2x2<Precision>(d(ll), e(ll), d(ll+1), shift, r);
    70327043                }
    7033                
     7044
    70347045                //
    70357046                // Test if shift negligible, and if so set to zero
    … …  
    70437054                }
    70447055            }
    7045            
     7056
    70467057            //
    70477058            // Increment iteration count
    70487059            //
    70497060            iter = iter+m-ll;
    7050            
     7061
    70517062            //
    70527063            // If SHIFT = 0, do simplified QR iteration
    … …  
    70567067                if( idir==1 )
    70577068                {
    7058                    
     7069
    70597070                    //
    70607071                    // Chase bulge from top to bottom
    … …  
    70807091                    d(m) = h*oldcs;
    70817092                    e(m-1) = h*oldsn;
    7082                    
     7093
    70837094                    //
    70847095                    // Update singular vectors
    … …  
    70967107                        rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work2, work3, c, ctemp);
    70977108                    }
    7098                    
     7109
    70997110                    //
    71007111                    // Test convergence
    … …  
    71077118                else
    71087119                {
    7109                    
     7120
    71107121                    //
    71117122                    // Chase bulge from bottom to top
    … …  
    71317142                    d(ll) = h*oldcs;
    71327143                    e(ll) = h*oldsn;
    7133                    
     7144
    71347145                    //
    71357146                    // Update singular vectors
    … …  
    71477158                        rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work0, work1, c, ctemp);
    71487159                    }
    7149                    
     7160
    71507161                    //
    71517162                    // Test convergence
    … …  
    71597170            else
    71607171            {
    7161                
     7172
    71627173                //
    71637174                // Use nonzero shift
    … …  
    71657176                if( idir==1 )
    71667177                {
    7167                    
     7178
    71687179                    //
    71697180                    // Chase bulge from top to bottom
    … …  
    71987209                    }
    71997210                    e(m-1) = f;
    7200                    
     7211
    72017212                    //
    72027213                    // Update singular vectors
    … …  
    72147225                        rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work2, work3, c, ctemp);
    72157226                    }
    7216                    
     7227
    72177228                    //
    72187229                    // Test convergence
    … …  
    72257236                else
    72267237                {
    7227                    
     7238
    72287239                    //
    72297240                    // Chase bulge from bottom to top
    … …  
    72587269                    }
    72597270                    e(ll) = f;
    7260                    
     7271
    72617272                    //
    72627273                    // Test convergence
    … …  
    72667277                        e(ll) = 0;
    72677278                    }
    7268                    
     7279
    72697280                    //
    72707281                    // Update singular vectors if desired
    … …  
    72847295                }
    72857296            }
    7286            
     7297
    72877298            //
    72887299            // QR iteration finished, go back and check convergence
    … …  
    72907301            continue;
    72917302        }
    7292        
     7303
    72937304        //
    72947305        // All singular values converged, so make them positive
    … …  
    72997310            {
    73007311                d(i) = -d(i);
    7301                
     7312
    73027313                //
    73037314                // Change sign of singular vectors, if desired
    … …  
    73097320            }
    73107321        }
    7311        
     7322
    73127323        //
    73137324        // Sort the singular values into decreasing order (insertion sort on
    … …  
    73167327        for(i=1; i<=n-1; i++)
    73177328        {
    7318            
     7329
    73197330            //
    73207331            // Scan for smallest D(I)
    … …  
    73327343            if( isub!=n+1-i )
    73337344            {
    7334                
     7345
    73357346                //
    73367347                // Swap singular values and vectors
    … …  
    74357446                if( au==0 )
    74367447                {
    7437                    
     7448
    74387449                    //
    74397450                    // Avoid possible harmful underflow if exponent range
    … …  
    75007511        ht = h;
    75017512        ha = amp::abs<Precision>(h);
    7502        
     7513
    75037514        //
    75047515        // PMAX points to the maximum absolute element of matrix
    … …  
    75117522        if( swp )
    75127523        {
    7513            
     7524
    75147525            //
    75157526            // Now FA .ge. HA
    … …  
    75277538        if( ga==0 )
    75287539        {
    7529            
     7540
    75307541            //
    75317542            // Diagonal matrix
    … …  
    75467557                if( fa/ga<amp::ampf<Precision>::getAlgoPascalEpsilon() )
    75477558                {
    7548                    
     7559
    75497560                    //
    75507561                    // Case of very large GA
    … …  
    75707581            if( gasmal )
    75717582            {
    7572                
     7583
    75737584                //
    75747585                // Normal case
    … …  
    76017612                if( mm==0 )
    76027613                {
    7603                    
     7614
    76047615                    //
    76057616                    // Note that M is very tiny
    … …  
    76407651            snr = srt;
    76417652        }
    7642        
     7653
    76437654        //
    76447655        // Correct signs of SSMAX and SSMIN
Note: See TracChangeset for help on using the changeset viewer.