Wrong with my codes according to sgemv

Open discussion for MAGMA

Wrong with my codes according to sgemv

Postby cz0397 » Wed Apr 10, 2013 4:09 am

Hi.I attampt to compile Jacobi iteration by magma. However, I meet some problems.First, I want to know whether magma has given some method to solve Ax = b.
Next, according to the testing in magma, I write codes to solve Ax = b. The result is not correct, and then I use gdb to debug my codes. I find such a problem that some codes are passed without being run.

Here is the total codes.

Code: Select all
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <cuda_runtime_api.h>
#include <cublas.h>

#include "flops.h"
#include "magma.h"
#include "magma_lapack.h"
#include "testings.h"
int main(int argc, char **argv)
{       
    TESTING_CUDA_INIT();
    magma_int_t row = 100;
    magma_int_t col = 100;
    magma_int_t szeA, szeb, szex;
    magma_int_t ione     = 1;
    magma_int_t ISEED[4] = {0,0,0,1};
    float *A, *M, *x, *b, *g;
    magma_int_t lda = ((row+31)/32)*32; //times of 32
    magma_int_t incx = 1;
    magma_int_t incb = 1;
   
    szeA = lda*col;
    szeb = row;
    szex = row;
    TESTING_MALLOC( A, float, szeA );
    TESTING_MALLOC( M, float, szeA );
    TESTING_MALLOC( x, float, szex );
    TESTING_MALLOC( b, float, szeb );
    TESTING_MALLOC( g, float, szeb );
    // generate rand matrix and vector

    //lapackf77_slarnv( &ione, ISEED, &szeA, A );
    //lapackf77_slarnv( &ione, ISEED, &szeb, b );
    for(int i = 0;i < row; ++i)
       for(int j = 0;j < col; ++j)
          A[i*col + j] = 0;
    for(int i = 0;i < row; ++i)
    {
       A[i*col + i] = 1;
       b[i] = 1;
    }

    // generate the M and g putting in A and b

    for(int i = 0;i < row; ++i)
       for(int j = 0;j < col; ++j)
          M[i*col + j] = A[i*col + j]/A[i*col + i];//D**(-1)*(L+U)
    for(int i = 0;i < row; ++i)
    {
       g[i] = b[i] / A[i*col + i];  //D**(-1)*b
       M[i*col + i] = 0;
       x[i] = 0;                    //initial vector
    }
    x[0] = 1;
    /*for(int i = 0;i < row; ++i)
    {
       for(int j = 0;j < col; ++j)
          printf("%f\t",A[i*col+j]);
       printf("\n");
    }*/
    //runing on GPU
    float *dA, *dM, *db, *dg, *dx;
    TESTING_DEVALLOC( dA, float, szeA );
    TESTING_DEVALLOC( dM, float, szeA );
    TESTING_DEVALLOC( db, float, szeb );
    TESTING_DEVALLOC( dg, float, szeb );
    TESTING_DEVALLOC( dx, float, szex );
    float error, work[1];
    magma_ssetmatrix( row, col, A, lda, dA, lda );
    magma_ssetmatrix( row, col, M, lda, dM, lda );

    char trans = MagmaNoTrans;
    float alpha = 1;
    float beta  = 1;
    magma_timestr_t  start, end;
    int num = 1;
    start = get_current_time();
    while(error < 1e-6 && num < 101)
    {
         magma_ssetvector( row, g, incb, dg, incb );//recover dg
         magma_sgetvector( row, dg, incb, g, incb );
         magmablas_sgemv( trans, row, col, alpha, dM, lda, dx, incx, beta, dg, incb );//the result(Mx+g) return in dg on GPU
         magma_sgetvector( row, dg, incb, x, incb );//put the result in x on CPU
         magma_ssetvector( col, x, incx, dx, incx );//update dx
         alpha = -1;
         magma_ssetvector( row, b, incb, db, incb );//recover db
         magmablas_sgemv( trans, row, col, alpha, dA, lda, dx, incx, beta, db, incb );//the result(b-Ax) return in db on GPU
         magma_sgetvector( row, db, incb, x, incb );//put the result in x on CPU
         error = lapackf77_slange( "M", &row, &ione, x, &row, work );
         num++;     
    }
    end = get_current_time();
    float time = GetTimerValue(start, end);
    printf("Error is %f\n",error);
    printf("Time is %f\n",time);
    /* Free Memory */
    TESTING_FREE( A );
    TESTING_FREE( M );
    TESTING_FREE( x );
    TESTING_FREE( b );
    TESTING_FREE( g );

    TESTING_DEVFREE( dA );
    TESTING_DEVFREE( dM );
    TESTING_DEVFREE( dx );
    TESTING_DEVFREE( db );
    TESTING_DEVFREE( dg );

    /* Free device */
    TESTING_CUDA_FINALIZE();
    return EXIT_SUCCESS;
}


The following meets the problem.

Code: Select all
while(error < 1e-6 && num < 101)
    {
         magma_ssetvector( row, g, incb, dg, incb );//recover dg
         magma_sgetvector( row, dg, incb, g, incb );
         magmablas_sgemv( trans, row, col, alpha, dM, lda, dx, incx, beta, dg, incb );//the result(Mx+g) return in dg on GPU
         magma_sgetvector( row, dg, incb, x, incb );//put the result in x on CPU
         magma_ssetvector( col, x, incx, dx, incx );//update dx
         alpha = -1;
         magma_ssetvector( row, b, incb, db, incb );//recover db
         magmablas_sgemv( trans, row, col, alpha, dA, lda, dx, incx, beta, db, incb );//the result(b-Ax) return in db on GPU
         magma_sgetvector( row, db, incb, x, incb );//put the result in x on CPU
         error = lapackf77_slange( "M", &row, &ione, x, &row, work );
         num++;     
    }


I always do the 3rd row, then going to the last. and it will only do two loop.
Code: Select all
gdb ./jacobi1 GNU gdb (Ubuntu/Linaro 7.3-0ubuntu2) 7.3-2011.08
Copyright (C) 2011 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "i686-linux-gnu".
For bug reporting instructions, please see:
<http://bugs.launchpad.net/gdb-linaro/>...
Reading symbols from /home/mathlab/magma-1.3.0/testing/jacobi1...done.
(gdb) b 83
Breakpoint 1 at 0x80493b5: file jacobi1.cpp, line 83.
(gdb) r
Starting program: /home/mathlab/magma-1.3.0/testing/jacobi1
[Thread debugging using libthread_db enabled]
[New Thread 0x5d27b70 (LWP 3569)]
MAGMA 1.3.0
device 0: GeForce 405, 1402.0 MHz clock, 511.6 MB memory, capability 1.2

Breakpoint 1, main (argc=Cannot access memory at address 0x0
) at jacobi1.cpp:83
83            magma_ssetvector( row, g, incb, dg, incb );//recover dg
(gdb) n
93            num++;     
(gdb) p num
$1 = <optimized out>
(gdb) n
83            magma_ssetvector( row, g, incb, dg, incb );//recover dg
(gdb) p num
$2 = <optimized out>
(gdb) n
84            magma_sgetvector( row, dg, incb, g, incb );
(gdb) p g[0]
$3 = 1
(gdb) p g[1]
$4 = 1
(gdb) n
85            magmablas_sgemv( trans, row, col, alpha, dM, lda, dx, incx, beta, dg, incb );//the result(Mx+g) return in dg on GPU
(gdb) p g[0]
$5 = 1
(gdb) p g[1]
$6 = 1
(gdb) p g[99]
$7 = 1
(gdb) n
90            magmablas_sgemv( trans, row, col, alpha, dA, lda, dx, incx, beta, db, incb );//the result(b-Ax) return in db on GPU
(gdb) n
85            magmablas_sgemv( trans, row, col, alpha, dM, lda, dx, incx, beta, dg, incb );//the result(Mx+g) return in dg on GPU
(gdb) p num
$8 = <optimized out>
(gdb) p error
$9 = <optimized out>
(gdb) n
86            magma_sgetvector( row, dg, incb, x, incb );//put the result in x on CPU
(gdb) n
87            magma_ssetvector( col, x, incx, dx, incx );//update dx
(gdb) n
89            magma_ssetvector( row, b, incb, db, incb );//recover db
(gdb) n
90            magmablas_sgemv( trans, row, col, alpha, dA, lda, dx, incx, beta, db, incb );//the result(b-Ax) return in db on GPU
(gdb) n
91            magma_sgetvector( row, db, incb, x, incb );//put the result in x on CPU
(gdb) n
92            error = lapackf77_slange( "M", &row, &ione, x, &row, work );
(gdb) n
81       while(error < 1e-6 && num < 101)
(gdb) n
95       end = get_current_time();
(gdb) p error
$10 = <optimized out>
(gdb) p num
$11 = 2
(gdb) n
96       float time = GetTimerValue(start, end);
(gdb) n
97       printf("Error is %f\n",error);
(gdb) n
Error is inf


I think there must be some wrongs with my codes. Please help me to find it, Thanks a lot!
cz0397
 
Posts: 15
Joined: Tue Feb 26, 2013 10:22 pm

Re: Wrong with my codes according to sgemv

Postby mgates3 » Wed Apr 10, 2013 6:15 pm

The routines magma_[sdcz]gesv will solve Ax=b for a dense, nonsymmetric matrix A.

I didn't follow exactly what you expected your code to do versus what it actually did. The best I can do is offer some observations:

- I don't understand why you do setvector and then immediately getvector. The getvector won't change anything. That is, it seems you are doing:
dg = g;
g = dg; // not needed
Also later doing getvector followed by setvector.

- You change alpha in the loop, but don't change it back. Better to use descriptive constants like c_zero, c_one, c_neg_one.

- Your loop termination looks wrong. Shouldn't error < 1e-6 be > 1e-6? Also, you need to initialize error before using it in the loop condition. Using gcc -O1 -Wuninitialized can sometimes find this problem.

-mark
mgates3
 
Posts: 427
Joined: Fri Jan 06, 2012 2:13 pm

Re: Wrong with my codes according to sgemv

Postby cz0397 » Fri Apr 12, 2013 12:34 am

Hi, mark.

First, I decide to solve Ax = b, and use Jacobi iteration, which is x = M*x + g. You must already know. And you know the sgemv does y = alpha*A*x + beta*y, and then this will change my g. Thus I have to do getvector followed by setvector, which is the second. The first point in my code is just to test whether I transfer correctly. The problem is that I transfer data too often.

Code: Select all
81       start = get_current_time();
(gdb) n
80       int num = 1;
(gdb) n
93            error = lapackf77_slange( "M", &row, &ione, x, &row, work );
(gdb) p error
$10 = 1
(gdb) n
85            magma_ssetvector( row, g, incb, dg, incb );//recover dg
(gdb) p error
$11 = <optimized out>
(gdb) p x[0]
$12 = 1
(gdb) p x[1]
$13 = 0
(gdb) p x[2]
$14 = 0
(gdb) p x[20]
$15 = 0
(gdb) p x[99]
$16 = 0

The program will do the 93th first, which change value of error. But the norm of x should be one. Do I use lapackf77_slange wrongly?
And I meet similar problem which I get the value from GPU.
Code: Select all
87            magma_sgetvector( row, dg, incb, x, incb );//put the result in x on CPU
(gdb) p x[0]
$18 = 1
(gdb) n
88            magma_ssetvector( col, x, incx, dx, incx );//update dx
(gdb) p x[1]
$19 = nan(0x7fffff)



Second, the alpha should be change back, and error should be given an initial value.
Shouldn't error < 1e-6 be > 1e-6?
I forget the 'while' I use.


Thanks for helping!
cz0397
 
Posts: 15
Joined: Tue Feb 26, 2013 10:22 pm

Re: Wrong with my codes according to sgemv

Postby cz0397 » Sun Apr 14, 2013 11:47 pm

Hi mark.

I get the problem as said in last post.
the norm of x should be one. Do I use lapackf77_slange wrongly?
And I meet similar problem which I get the value from GPU. Next I run the test file, and find that all function goes correctly. I know there must be some problems in my codes, but I have no idea where they are.
Code: Select all
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <cuda_runtime_api.h>
#include <cublas.h>

#include "flops.h"
#include "magma.h"
#include "magma_lapack.h"
#include "testings.h"
int main(int argc, char **argv)
{       
    TESTING_CUDA_INIT();
    magma_int_t row = 100;
    magma_int_t col = 100;
    magma_int_t szeA, szeb, szex;
    magma_int_t ione     = 1;
    magma_int_t ISEED[4] = {0,0,0,1};
    float *A, *M, *x, *b, *g;
    magma_int_t lda = ((row+31)/32)*32; //times of 32
    magma_int_t incx = 1;
    magma_int_t incb = 1;
   
    szeA = lda*col;
    szeb = row;
    szex = row;
    TESTING_MALLOC( A, float, szeA );
    TESTING_MALLOC( M, float, szeA );
    TESTING_MALLOC( x, float, szex );
    TESTING_MALLOC( b, float, szeb );
    TESTING_MALLOC( g, float, szeb );
    // generate rand matrix and vector

    //lapackf77_slarnv( &ione, ISEED, &szeA, A );
    //lapackf77_slarnv( &ione, ISEED, &szeb, b );
    for(int i = 0;i < row; ++i)
       for(int j = 0;j < col; ++j)
          A[i*col + j] = 0;
    for(int i = 0;i < row; ++i)
    {
       A[i*col + i] = 1;
       b[i] = 1;
    }

    // generate the M and g putting in A and b

    for(int i = 0;i < row; ++i)
       for(int j = 0;j < col; ++j)
          M[i*col + j] = A[i*col + j]/A[i*col + i];//D**(-1)*(L+U)
    for(int i = 0;i < row; ++i)
    {
       g[i] = b[i] / A[i*col + i];  //D**(-1)*b
       M[i*col + i] = 0;
       x[i] = 0;                    //initial vector
    }
    x[0] = 1;
    /*for(int i = 0;i < row; ++i)
    {
       for(int j = 0;j < col; ++j)
          printf("%f\t",A[i*col+j]);
       printf("\n");
    }*/
    //runing on GPU
    float *dA, *dM, *db, *dg, *dx;
    TESTING_DEVALLOC( dA, float, szeA );
    TESTING_DEVALLOC( dM, float, szeA );
    TESTING_DEVALLOC( db, float, szeb );
    TESTING_DEVALLOC( dg, float, szeb );
    TESTING_DEVALLOC( dx, float, szex );
    float error = MAGMA_S_ONE;
    float work[1];
    magma_ssetmatrix( row, col, A, lda, dA, lda );
    magma_ssetmatrix( row, col, M, lda, dM, lda );

    char trans = MagmaNoTrans;
    float alpha;
    float beta  = MAGMA_S_ONE;
    magma_timestr_t  start, end;
    int num = 1;
    start = get_current_time();
    while(error > 1e-6 && num < 101)
    {
         //solve x
         alpha = MAGMA_S_ONE;
         magma_ssetvector( row, g, incb, dg, incb );//recover dg
         magmablas_sgemv( trans, row, col, alpha, dM, lda, dx, incx, beta, dg, incb );//the result(Mx+g) return in dg on GPU
         magma_sgetvector( row, dg, incb, x, incb );//put the result in x on CPU
         //solve error
         magma_ssetvector( col, x, incx, dx, incx );//update dx
         alpha = MAGMA_S_NEG_ONE;
         magma_ssetvector( row, b, incb, db, incb );//recover db
         magmablas_sgemv( trans, row, col, alpha, dA, lda, dx, incx, beta, db, incb );//the result(b-Ax) return in db on GPU
         magma_sgetvector( row, db, incb, x, incb );//put the result in x on CPU
         error = lapackf77_slange( "M", &row, &ione, x, &row, work );
         num++;     
    }
    end = get_current_time();
    float time = GetTimerValue(start, end);
    printf("Error is %f\n",error);
    printf("Time is %f\n",time);
    /* Free Memory */
    TESTING_FREE( A );
    TESTING_FREE( M );
    TESTING_FREE( x );
    TESTING_FREE( b );
    TESTING_FREE( g );

    TESTING_DEVFREE( dA );
    TESTING_DEVFREE( dM );
    TESTING_DEVFREE( dx );
    TESTING_DEVFREE( db );
    TESTING_DEVFREE( dg );

    /* Free device */
    TESTING_CUDA_FINALIZE();
    return EXIT_SUCCESS;
}

Thank you!
cz0397
 
Posts: 15
Joined: Tue Feb 26, 2013 10:22 pm


Return to User discussion

Who is online

Users browsing this forum: Bing [Bot] and 1 guest

cron