matrix multiplication on parallella 16 core.

Any technical questions about the Epiphany chip and Parallella HW Platform.

Moderator: aolofsson

matrix multiplication on parallella 16 core.

Postby sriranjan » Wed Feb 03, 2016 7:49 am

I want to buy a parallella board. I want to run my program which involves lot of matrix multiplication.How do I do it like yaniv showed in the video?
https://www.youtube.com/watch?v=DkctH7_tYSc
Could someone point me to the exact library and example code to do matrix multiplication on host and the parallella board. Right now my algo takes 120 ms and if I run all the matrix multiplications on the t16 core board, how fast will my program execute? Please let me know.
sriranjan
 
Posts: 12
Joined: Mon Feb 01, 2016 3:17 pm

Re: matrix multiplication on parallella 16 core.

Postby ubii » Wed Feb 03, 2016 9:25 pm

User avatar
ubii
 
Posts: 71
Joined: Sun Dec 16, 2012 7:18 pm
Location: US

Re: matrix multiplication on parallella 16 core.

Postby sriranjan » Thu Feb 04, 2016 10:30 am

Thanks I will try this.
sriranjan
 
Posts: 12
Joined: Mon Feb 01, 2016 3:17 pm

Re: matrix multiplication on parallella 16 core.

Postby sriranjan » Thu Feb 04, 2016 10:36 am

Where is yaniv's example code. That example was like just one page. Could someone point me to yaniv's matrix multiplication code?
sriranjan
 
Posts: 12
Joined: Mon Feb 01, 2016 3:17 pm

Re: matrix multiplication on parallella 16 core.

Postby sriranjan » Thu Feb 04, 2016 10:42 am

My code is filled with matrix computations of size 4000x2. I am posting a function that I use in my code below.Can this be parallelized or how shall I make this execute on the coprocessor. Someone suggested me to use MPI with ubuntu 14.04 which would parallelize the code by itself.
Code: Select all
void svdcmp(mat A, int M, int N, vect W, mat V)
{
    /* Householder reduction to bidiagonal form. */
    int NM;
    scal C;
    scal F;
    scal G = 0.0;
    scal H;
    scal S;
    scal X;
    scal Y;
    scal Z;
    scal Scale = 0.0;
    scal ANorm = 0.0;
    scal tmp;
    int flag;
    int i;
    int its;
    int j;
    int jj;
    int k;
    int l;
   vect rv1;


   rv1 = vect_create(N);

    if( M < N ) {
        fprintf( stderr, "You must augment A with extra zero rows.\n" );
        return;
    }

    for( i = 0; i < N; ++i ) {
        l = i + 1;
        rv1[i] = Scale * G;
        G = 0.0;
        S = 0.0;
        Scale = 0.0;
        if( i < M ) {
            for( k = i; k < M; ++k ) {
                Scale = Scale + fabs( A[k][i] );
            }
            if( Scale != 0.0 ) {
                for( k = i; k < M; ++k ) {
                    A[k][i] = A[k][i] / Scale;
                    S = S + A[k][i] * A[k][i];
                }
                F = A[i][i];
                G = sqrt(S);
                if( F > 0.0 ) {
                    G = -G;
                }
                H = F * G - S;
                A[i][i] = F - G;
                if( i != (N-1) ) {
                    for( j = l; j < N; ++j ) {
                        S = 0.0;
                        for( k = i; k < M; ++k ) {
                            S = S + A[k][i] * A[k][j];
                        }
                        F = S / H;
                        for( k = i; k < M; ++k ) {
                            A[k][j] = A[k][j] + F * A[k][i];
                        }
                    }
                }
                for( k = i; k < M; ++k ) {
                    A[k][i] = Scale * A[k][i];
                }
            }
        }

        W[i] = Scale * G;
        G = 0.0;
        S = 0.0;
        Scale = 0.0;
        if( (i < M) && (i != (N-1)) ) {
            for( k = l; k < N; ++k ) {
                Scale = Scale + fabs( A[i][k] );
            }
            if( Scale != 0.0 ) {
                for( k = l; k < N; ++k ) {
                    A[i][k] = A[i][k] / Scale;
                    S = S + A[i][k] * A[i][k];
                }
                F = A[i][l];
                G = sqrt(S);
                if( F > 0.0 ) {
                    G = -G;
                }
                H = F * G - S;
                A[i][l] = F - G;
                for( k = l; k < N; ++k ) {
                    rv1[k] = A[i][k] / H;
                }
                if( i != (M-1) ) {
                    for( j = l; j < M; ++j ) {
                        S = 0.0;
                        for( k = l; k < N; ++k ) {
                            S = S + A[j][k] * A[i][k];
                        }
                        for( k = l; k < N; ++k ) {
                            A[j][k] = A[j][k] + S * rv1[k];
                        }
                    }
                }
                for( k = l; k < N; ++k ) {
                    A[i][k] = Scale * A[i][k];
                }
            }
        }
        tmp = fabs( W[i] ) + fabs( rv1[i] );
        if( tmp > ANorm )
            ANorm = tmp;
    }

    /* Accumulation of right-hand transformations. */
    for( i = N-1; i >= 0; --i ) {
        if( i < (N-1) ) {
            if( G != 0.0 ) {
                for( j = l; j < N; ++j ) {
                    V[j][i] = (A[i][j] / A[i][l]) / G;
                }
                for( j = l; j < N; ++j ) {
                    S = 0.0;
                    for( k = l; k < N; ++k ) {
                        S = S + A[i][k] * V[k][j];
                    }
                    for( k = l; k < N; ++k ) {
                        V[k][j] = V[k][j] + S * V[k][i];
                    }
                }
            }
            for( j = l; j < N; ++j ) {
                V[i][j] = 0.0;
                V[j][i] = 0.0;
            }
        }
        V[i][i] = 1.0;
        G = rv1[i];
        l = i;
    }

    /* Accumulation of left-hand transformations. */
    for( i = N-1; i >= 0; --i ) {
        l = i + 1;
        G = W[i];
        if( i < (N-1) ) {
            for( j = l; j < N; ++j ) {
                A[i][j] = 0.0;
            }
        }
        if( G != 0.0 ) {
            G = 1.0 / G;
            if( i != (N-1) ) {
                for( j = l; j < N; ++j ) {
                    S = 0.0;
                    for( k = l; k < M; ++k ) {
                        S = S + A[k][i] * A[k][j];
                    }
                    F = (S / A[i][i]) * G;
                    for( k = i; k < M; ++k ) {
                        A[k][j] = A[k][j] + F * A[k][i];
                    }
                }
            }
            for( j = i; j < M; ++j ) {
                A[j][i] = A[j][i] * G;
            }
        } else {
            for( j = i; j < M; ++j ) {
                A[j][i] = 0.0;
            }
        }
        A[i][i] = A[i][i] + 1.0;
    }

    /* Diagonalization of the bidiagonal form.
       Loop over singular values. */
    for( k = (N-1); k >= 0; --k ) {
        /* Loop over allowed iterations. */
        for( its = 1; its <= 30; ++its ) {
            /* Test for splitting.
               Note that rv1[0] is always zero. */
            flag = true;
            for( l = k; l >= 0; --l ) {
                NM = l - 1;
                if( (fabs(rv1[l]) + ANorm) == ANorm ) {
                    flag = false;
                    break;
                } else if( (fabs(W[NM]) + ANorm) == ANorm ) {
                    break;
                }
            }

            /* Cancellation of rv1[l], if l > 0; */
            if( flag ) {
                C = 0.0;
                S = 1.0;
                for( i = l; i <= k; ++i ) {
                    F = S * rv1[i];
                    if( (fabs(F) + ANorm) != ANorm ) {
                        G = W[i];
                        H = sqrt( F * F + G * G );
                        W[i] = H;
                        H = 1.0 / H;
                        C = ( G * H );
                        S = -( F * H );
                        for( j = 0; j < M; ++j ) {
                            Y = A[j][NM];
                            Z = A[j][i];
                            A[j][NM] = (Y * C) + (Z * S);
                            A[j][i] = -(Y * S) + (Z * C);
                        }
                    }
                }
            }
            Z = W[k];
            /* Convergence. */
            if( l == k ) {
                /* Singular value is made nonnegative. */
                if( Z < 0.0 ) {
                    W[k] = -Z;
                    for( j = 0; j < N; ++j ) {
                        V[j][k] = -V[j][k];
                    }
                }
                break;
            }

            if( its >= 30 ) {
                fprintf( stderr, "No convergence in 30 iterations.\n" );
                return;
            }

            X = W[l];
            NM = k - 1;
            Y = W[NM];
            G = rv1[NM];
            H = rv1[k];
            F = ((Y-Z)*(Y+Z) + (G-H)*(G+H)) / (2.0*H*Y);
            G = sqrt( F * F + 1.0 );
            tmp = G;
            if( F < 0.0 )
                tmp = -tmp;
            F = ((X-Z)*(X+Z) + H*((Y/(F+tmp))-H)) / X;

            /* Next QR transformation. */
            C = 1.0;
            S = 1.0;
            for( j = l; j <= NM; ++j ) {
                i = j + 1;
                G = rv1[i];
                Y = W[i];
                H = S * G;
                G = C * G;
                Z = sqrt( F * F + H * H );
                rv1[j] = Z;
                C = F / Z;
                S = H / Z;
                F = (X * C) + (G * S);
                G = -(X * S) + (G * C);
                H = Y * S;
                Y = Y * C;
                for( jj = 0; jj < N; ++jj ) {
                    X = V[jj][j];
                    Z = V[jj][i];
                    V[jj][j] = (X * C) + (Z * S);
                    V[jj][i] = -(X * S) + (Z * C);
                }
                Z = sqrt( F * F + H * H );
                W[j] = Z;

                /* Rotation can be arbitrary if Z = 0. */
                if( Z != 0.0 ) {
                    Z = 1.0 / Z;
                    C = F * Z;
                    S = H * Z;
                }
                F = (C * G) + (S * Y);
                X = -(S * G) + (C * Y);
                for( jj = 0; jj < M; ++jj ) {
                    Y = A[jj][j];
                    Z = A[jj][i];
                    A[jj][j] = (Y * C) + (Z * S);
                    A[jj][i] = -(Y * S) + (Z * C);
                }
            }
            rv1[l] = 0.0;
            rv1[k] = F;
            W[k] = X;
        }
    }
   vect_delete(rv1);
   
    return;
}

sriranjan
 
Posts: 12
Joined: Mon Feb 01, 2016 3:17 pm

Re: matrix multiplication on parallella 16 core.

Postby sriranjan » Thu Feb 04, 2016 1:54 pm

How do I convert my code to work on the 16 core. Shall I download OpenMP and it might parallelize it?
sriranjan
 
Posts: 12
Joined: Mon Feb 01, 2016 3:17 pm

Re: matrix multiplication on parallella 16 core.

Postby sriranjan » Thu Feb 04, 2016 2:49 pm

The following are the functions that I need to run on the 16 core chip. How shall I parallelize them?
Code: Select all
void mat_delete(mat M, int rows, int cols);
vect vect_create(int n);
void vect_delete(vect v);
void mat_zeroize(mat M, int rows, int cols);
void vect_zeroize(vect v, int n);

void mat_copy(mat M, int rows, int cols, mat Md);
void vect_apply_fx(vect v, int n, scal (*fx)(), scal par);
void mat_apply_fx(mat M, int rows, int cols, scal (*fx)(), scal par);
void mat_mean_rows(mat M, int rows, int cols, vect v);
scal mat_max_diag(mat M, int rows, int cols);
void mat_diag(vect v, int n, mat R);
void mat_transpose(mat M, int rows, int cols, mat R);
void mat_inverse(mat M, int dim, mat R);
void mat_sub(mat A, mat B, int rows, int cols, mat R);
void mat_mult(mat A, int rows_A, int cols_A, mat B, int rows_B, int cols_B, mat R);
void mat_center(mat M, int rows, int cols, vect means);
void mat_decenter(mat M, int rows, int cols, vect means);
sriranjan
 
Posts: 12
Joined: Mon Feb 01, 2016 3:17 pm

Re: matrix multiplication on parallella 16 core.

Postby ubii » Thu Feb 04, 2016 6:05 pm

I was not able to find the original matmul-16 on github, but listed below is the original white paper written by Yaniv, with a link to the source code.

http://www.adapteva.com/white-papers/sc ... -matrices/
User avatar
ubii
 
Posts: 71
Joined: Sun Dec 16, 2012 7:18 pm
Location: US

Re: matrix multiplication on parallella 16 core.

Postby sriranjan » Fri Feb 05, 2016 2:44 am

So how do I port the rest of my code to fit on the parallella?
sriranjan
 
Posts: 12
Joined: Mon Feb 01, 2016 3:17 pm


Return to Epiphany and Parallella Q & A

Who is online

Users browsing this forum: No registered users and 15 guests

cron