Function rmatrixgemm

src/ablas.cpp:1521–1615

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

     Bochkanov Sergey
*************************************************************************/
void rmatrixgemm(int m,
     int n,
     int k,
     double alpha,
     const ap::real_2d_array& a,
     int ia,
     int ja,
     int optypea,
     const ap::real_2d_array& b,
     int ib,
     int jb,
     int optypeb,
     double beta,
     ap::real_2d_array& c,
     int ic,
     int jc)
{
    int s1;
    int s2;
    int bs;

    bs = ablasblocksize(a);
    if( m<=bs&&n<=bs&&k<=bs )
    {
        rmatrixgemmk(m, n, k, alpha, a, ia, ja, optypea, b, ib, jb, optypeb, beta, c, ic, jc);
        return;
    }
    if( m>=n&&m>=k )
    {
        
        //
        // A*B = (A1 A2)^T*B
        //
        ablassplitlength(a, m, s1, s2);
        if( optypea==0 )
        {
            rmatrixgemm(s1, n, k, alpha, a, ia, ja, optypea, b, ib, jb, optypeb, beta, c, ic, jc);
            rmatrixgemm(s2, n, k, alpha, a, ia+s1, ja, optypea, b, ib, jb, optypeb, beta, c, ic+s1, jc);
        }
        else
        {
            rmatrixgemm(s1, n, k, alpha, a, ia, ja, optypea, b, ib, jb, optypeb, beta, c, ic, jc);
            rmatrixgemm(s2, n, k, alpha, a, ia, ja+s1, optypea, b, ib, jb, optypeb, beta, c, ic+s1, jc);
        }
        return;
    }
    if( n>=m&&n>=k )
    {
        
        //
        // A*B = A*(B1 B2)
        //
        ablassplitlength(a, n, s1, s2);
        if( optypeb==0 )
        {
            rmatrixgemm(m, s1, k, alpha, a, ia, ja, optypea, b, ib, jb, optypeb, beta, c, ic, jc);
            rmatrixgemm(m, s2, k, alpha, a, ia, ja, optypea, b, ib, jb+s1, optypeb, beta, c, ic, jc+s1);
        }