Function rmatrixsyrk

src/ablas.cpp:1295–1372

Same as CMatrixSYRK, but for real matrices OpType may be only 0 or 1. -- ALGLIB routine -- 16.12.2009 Bochkanov Sergey *************************************************************************/

     Bochkanov Sergey
*************************************************************************/
void rmatrixsyrk(int n,
     int k,
     double alpha,
     const ap::real_2d_array& a,
     int ia,
     int ja,
     int optypea,
     double beta,
     ap::real_2d_array& c,
     int ic,
     int jc,
     bool isupper)
{
    int s1;
    int s2;
    int bs;

    bs = ablasblocksize(a);
    if( n<=bs&&k<=bs )
    {
        rmatrixsyrk2(n, k, alpha, a, ia, ja, optypea, beta, c, ic, jc, isupper);
        return;
    }
    if( k>=n )
    {
        
        //
        // Split K
        //
        ablassplitlength(a, k, s1, s2);
        if( optypea==0 )
        {
            rmatrixsyrk(n, s1, alpha, a, ia, ja, optypea, beta, c, ic, jc, isupper);
            rmatrixsyrk(n, s2, alpha, a, ia, ja+s1, optypea, 1.0, c, ic, jc, isupper);
        }
        else
        {
            rmatrixsyrk(n, s1, alpha, a, ia, ja, optypea, beta, c, ic, jc, isupper);
            rmatrixsyrk(n, s2, alpha, a, ia+s1, ja, optypea, 1.0, c, ic, jc, isupper);
        }
    }
    else
    {
        
        //
        // Split N
        //
        ablassplitlength(a, n, s1, s2);
        if( optypea==0&&isupper )
        {
            rmatrixsyrk(s1, k, alpha, a, ia, ja, optypea, beta, c, ic, jc, isupper);
            rmatrixgemm(s1, s2, k, alpha, a, ia, ja, 0, a, ia+s1, ja, 1, beta, c, ic, jc+s1);
            rmatrixsyrk(s2, k, alpha, a, ia+s1, ja, optypea, beta, c, ic+s1, jc+s1, isupper);
            return;
        }
        if( optypea==0&&!isupper )
        {
            rmatrixsyrk(s1, k, alpha, a, ia, ja, optypea, beta, c, ic, jc, isupper);