Complex bicgstab on gpu

305 views
Skip to first unread message

Elliot

unread,
Nov 8, 2010, 8:55:01 PM11/8/10
to cusp-users
Hi,

I'm new to gpu programming and am trying to use cusp to speed up some
large complex matrix solves for matrices in csr format. As a first
step I had encouraging results in double precision for the poisson
equation example. For a 1000 x 1000 grid I was seeing a factor of 8
speedup on the gpu. Now I'm trying to get things working on a simple
complex example. I'm not 100% sure what has been implemented for
complex.

For my simple example I'm trying to solve the
coordinate_complex_general matrix (\cusp-library\testing\data\test)
example using a right hand side of 1+i. When I try running it on the
host I get a result which agrees with Matlab. When I try running it on
the gpu the program crashes with both the diagonal and identity
preconditioners. I'm compiling it with -arch sm_13.

Any suggestions would be greatly appreciated

Thanks

Elliot

#include <cusp/precond/diagonal.h>
#include <cusp/dia_matrix.h>
#include <cusp/csr_matrix.h>
#include <cusp/io/matrix_market.h>
#include <cusp/krylov/bicgstab.h>
#include <iostream>
#include <cusp/complex.h>
#include <cusp/blas.h>

typedef cusp::device_memory DM; // define Device Memory
typedef cusp::host_memory HM; // define Host Memory
typedef cusp::complex<double> Complex;

int main(void)
{
cusp::csr_matrix<int, Complex, HM> mycsr; // allocate memory
cusp::io::read_matrix_market_file(mycsr, "A.mtx"); //get a
csr_matrix in host memory

cusp::dia_matrix<int,Complex,HM>A(mycsr); //switch to dia_matrix
format
cusp::array1d<Complex, HM> x_host(A.num_rows, Complex(0,0));
cusp::array1d<Complex, HM> b_host(A.num_rows, Complex(1,1));
cusp::verbose_monitor<Complex> monitor_host(b_host, 100, 1e-5);
cusp::precond::diagonal<Complex, HM> M_host(A);
cusp::krylov::bicgstab(A, x_host, b_host, monitor_host, M_host);
cusp::print_matrix(x_host);

cusp::dia_matrix<int,Complex,DM>D(mycsr); //have dia_matrix in
device memory
cusp::array1d<Complex, DM> x(D.num_rows, Complex(0,0));
cusp::array1d<Complex, DM> b(D.num_rows, Complex(1,1));
cusp::verbose_monitor<Complex> monitor(b, 100, 1e-5);
cusp::precond::diagonal<Complex, DM> M(D);
cusp::krylov::bicgstab(D, x, b, monitor, M);

return 0;
}

Filipe Maia

unread,
Nov 8, 2010, 10:59:08 PM11/8/10
to cusp-...@googlegroups.com
Hi,

Here the device part of the code doesn't even compile.
I get:

./cusp/detail/device/spmv/dia.h(87): error: ambiguous "?" operation: second operand of type "int" can be converted to third operand type "Complex", and vice versa

What CUDA/nvcc version do you have?

Cheers,
Filipe


--
You received this message because you are subscribed to the Google Groups "cusp-users" group.
To post to this group, send email to cusp-...@googlegroups.com.
To unsubscribe from this group, send email to cusp-users+...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/cusp-users?hl=en.


Filipe Maia

unread,
Nov 8, 2010, 11:04:42 PM11/8/10
to cusp-...@googlegroups.com
I opened an issue about it [1].
Try the patch I uploaded there.

Luca

unread,
Nov 9, 2010, 9:24:11 AM11/9/10
to cusp-users
I am also testing the complex version ... I really need it.

I can compile with matrices of type cusp::ell_matrix<int,
cusp::complex<float>, cusp::device_memory> and use the multiply
algorithm but the result is wrong. It looks like the imaginary part is
set to zero in the matrix. Is this related to the ongoing effort to
implement complex numbers?

Cheers,
Luca

ps: I use CUDA 3.2 with arch=sm_21

On Nov 9, 4:59 am, Filipe Maia <filipe.c.m...@gmail.com> wrote:
> Hi,
>
> Here the device part of the code doesn't even compile.
> I get:
>
> ./cusp/detail/device/spmv/dia.h(87): error: ambiguous "?" operation: second
> operand of type "int" can be converted to third operand type "Complex", and
> vice versa
>
> What CUDA/nvcc version do you have?
>
> Cheers,
> Filipe
>
> > cusp-users+...@googlegroups.com<cusp-users%2Bunsu...@googlegroups.com>
> > .

Luca

unread,
Nov 9, 2010, 11:30:23 AM11/9/10
to cusp-users
Sorry I forgot to mention that I first define a matrix
cusp::coo_matrix <int, cusp::complex<float>, cusp::host_memory> which
I then convert to ell_matrix on the device.

If I output the coo_matrix with cusp::io::write_matrix_market_file,
this is where the imaginary part is already lost.

Cheers,
Luca

Filipe Maia

unread,
Nov 9, 2010, 11:46:10 AM11/9/10
to cusp-...@googlegroups.com
On Tue, Nov 9, 2010 at 08:30, Luca <cov...@gmail.com> wrote:
Sorry I forgot to mention that I first define a matrix
cusp::coo_matrix <int, cusp::complex<float>, cusp::host_memory> which
I then convert to ell_matrix on the device.

If I output the coo_matrix with cusp::io::write_matrix_market_file,
this is where the imaginary part is already lost.


It's entirely possible that things don't work with the complex numbers as all the code was written assuming reals and the complex numbers were added later.
We need more people like you that try it  out and let the others know what things don't work so that they can be fixed.
 
To unsubscribe from this group, send email to cusp-users+...@googlegroups.com.

Filipe Maia

unread,
Nov 9, 2010, 2:51:13 PM11/9/10
to cusp-...@googlegroups.com
On Tue, Nov 9, 2010 at 08:30, Luca <cov...@gmail.com> wrote:
Sorry I forgot to mention that I first define a matrix
cusp::coo_matrix <int, cusp::complex<float>, cusp::host_memory> which
I then convert to ell_matrix on the device.

If I output the coo_matrix with cusp::io::write_matrix_market_file,
this is where the imaginary part is already lost.

Could you please post the code you are using?
 

To unsubscribe from this group, send email to cusp-users+...@googlegroups.com.

Elliot

unread,
Nov 9, 2010, 2:52:28 PM11/9/10
to cusp-users
Thanks for the quick reply and the patch. It seems like there is quite
a bit of interest in the complex implementation.

I downloaded your patch and tried re-running the code but it still
won't work.

The code will compile with some warnings. Two of them have to do with
dia.h

1>C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\bin/../
include\cusp/detail/device/spmv/dia.h(87): Warning: Cannot tell what
pointer points to, assuming global memory space
1>C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\bin/../
include\cusp/detail/device/spmv/dia.h(87): Warning: Cannot tell what
pointer points to, assuming global memory space

I'm trying to run this on a gts450 with I believe the newest version
of cuda

CUDA Driver Version: 3020
Device Number: 0
Device Name: GeForce GTS 450
Device Revision Number: 2.1
Global Memory Size: 1041694720

nvcc --version
nvcc: NVIDIA (R) Cuda compiler driver
Copyright (c) 2005-2010 NVIDIA Corporation
Built on Tue_Oct_19_02:27:10_PDT_2010
Cuda compilation tools, release 3.2, V0.2.1221

Cheers,

Elliot

Filipe Maia

unread,
Nov 9, 2010, 3:06:38 PM11/9/10
to cusp-...@googlegroups.com
On Tue, Nov 9, 2010 at 11:52, Elliot <elliot...@gmail.com> wrote:
Thanks for the quick reply and the patch. It seems like there is quite
a bit of interest in the complex implementation.

I downloaded your patch and tried re-running the code but it still
won't work.

The code will compile with some warnings. Two of them have to do with
dia.h

1>C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\bin/../
include\cusp/detail/device/spmv/dia.h(87): Warning: Cannot tell what
pointer points to, assuming global memory space
1>C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\bin/../
include\cusp/detail/device/spmv/dia.h(87): Warning: Cannot tell what
pointer points to, assuming global memory space

 
I can't reproduce these warnings on my machine. What compiler options did you use?
Also the warning seems related to the y[row] which was already there.

Elliot

unread,
Nov 9, 2010, 5:08:18 PM11/9/10
to cusp-users
So with your new patch the complex bicgstab example works for you now?

I do get a few warnings when I compile the code but nothing that
prevented linking. Things generally seemed to be working, and I didn't
understand the warnings, so I ignored the warnings ...
Here is the build and some of the warnings I get.

Cheers,

Elliot

1>Compiling with CUDA Build Rule...
1>"C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\\bin
\nvcc.exe" -arch sm_13 -ccbin "C:\Program Files (x86)\Microsoft
Visual Studio 8\VC\bin" -Xcompiler "/EHsc /W3 /nologo /O2 /Zi /
MT " -I"C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\
\include" -I"C:\ProgramData\NVIDIA Corporation\NVIDIA GPU Computing
SDK 3.2\C\common\inc" -maxrregcount=32 --compile -o "x64\Release
\sample.cu.obj" "c:\Users\Elliot\Documents\Visual Studio 2005\Projects
\CUDAWinApp2\CUDAWinApp2\sample.cu"
1>sample.cu

Some of the warnings I get

C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\bin/../include
\cusp/detail/device/spmv/dia.h(87): Warning: Cannot tell what pointer
points to, assuming global memory space

C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\bin/../include
\cusp/detail/device/spmv/dia.h(87): Warning: Cannot tell what pointer
points to, assuming global memory space

C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\bin/../include
\thrust/functional.h(409) : warning C4995: 'absolute_value': name was
marked as #pragma deprecated

C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\bin/../include
\cusp/detail/host/conversion.h(149) : warning C4267: 'argument' :
conversion from 'size_t' to 'int', possible loss of data

C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\bin/../include
\thrust/detail/dispatch/uninitialized_copy.h(46) : warning C4996:
'std::uninitialized_copy': Function call with parameters that may be
unsafe - this call relies on the caller to check that the passed
values are correct. To disable this warning, use -
D_SCL_SECURE_NO_WARNINGS. See documentation on how to use Visual C++
'Checked Iterators'

C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\bin/../include
\cusp/detail/host/conversion_utils.h(51) : warning C4267:
'initializing' : conversion from 'size_t' to 'IndexType', possible
loss of data

C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\bin/../include
\cusp/detail/host/conversion.h(514) : warning C4244: 'initializing' :
conversion from '__int64' to 'IndexType', possible loss of data

C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v3.2\bin/../include
\thrust/detail/dispatch/uninitialized_copy.h(46) : warning C4996:
'std::uninitialized_copy': Function call with parameters that may be
unsafe - this call relies on the caller to check that the passed
values are correct. To disable this warning, use -
D_SCL_SECURE_NO_WARNINGS. See documentation on how to use Visual C++
'Checked Iterators'
1> C:\Program Files (x86)\Microsoft Visual Studio 8\VC\include
\memory(99) : see declaration of 'std::uninitialized_copy'

caihon...@gmail.com

unread,
Jul 19, 2016, 2:09:55 AM7/19/16
to cusp-users, elliot...@gmail.com
Hello,

Can anyone send me a sample copy for the complex bicgstable code with CUSP? I write my own and it just doesn't work for the complex matrix.

Thanks.

Hongzhu

Steven Dalton

unread,
Jul 19, 2016, 10:00:08 AM7/19/16
to cusp-...@googlegroups.com
Hello,

  So you are able to compile your code but it fails at runtime? Can you send me a copy of your test matrix?

Steve

--
You received this message because you are subscribed to the Google Groups "cusp-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cusp-users+...@googlegroups.com.

To post to this group, send email to cusp-...@googlegroups.com.

caihon...@gmail.com

unread,
Jul 20, 2016, 1:31:47 AM7/20/16
to cusp-users
Hello Steve,

The attached is my code with two input mtx files for the matrix and right hand side. The code (complex number version) now works for bicgstab with diagonal preconditioner or without any preconditioner. However, I still get a lot of errors when I switch to the scaled_bridson_ainv preconditioner. You can try that commented line for the scaled_bridson_ainv preconditioner in my code. I have similar problems for other kinds of preconditioners. Could you please help me test it?

However, in order to use the diagonal preconditioner, I have to convert my A matrix from csr format to dia_matrix format. Otherwise, I get a lot of errors. I really don't know why.

In the meantime, how can I print the preconditioner matrix (M) use cusp::print? I can print other matrix. I get a lot of errors when I try to print M matrix in a similar way. 


I'm relative new to GPU and C++ Programming, my questions may looks stupid but I'm really looking forward to get helps from the community.

Thanks for your effort.

Regards.

Hongzhu
bicg_tmp_ReadMatrixRhs.cu
A_Matlab.mtx
b.mtx

Steven Dalton

unread,
Jul 22, 2016, 12:18:57 PM7/22/16
to cusp-...@googlegroups.com
Thanks for attaching your example code Hongzhu.

The scaled Bridson preconditioner needs to be updated to support complex numbers and tighter integration with thrust's design, it's on the list of issues [1].

I was able to compile your code using the smoothed_aggregation preconditioner. Although using complex types generates some warnings related to constructors in shared memory the code compiles fine, however SA does not currently support non-symmetric matrices.

The preconditioners do not support cusp::print since they are all a little different in the way they store and represent data. For SA you can use M.print() to display the relevant information about what the preconditioner has computed. 

You may receive some warnings about complex types but the diagonal preconditioner supports all sparse storage formats automatically. What version are you using?

Steve

caihon...@gmail.com

unread,
Jul 22, 2016, 1:26:48 PM7/22/16
to cusp-users
Hello Steve, 

Thanks for testing my code. As you mentioned that you can compile the code with smoothed_aggregation preconditioner. I have tried it again and it doesn't work in my machine actually. The attached is my error message.

Could you please send me the code that you modified to SA preconditioner and I will test it again. I guess I did the right thing to call SA preconditioner and it's really straightforward. I don't know what's the problem. 

I really appreciate for your help and I'm looking forward to hear from you soon.

Regards.

Hongzhu
message

caihon...@gmail.com

unread,
Jul 22, 2016, 2:00:04 PM7/22/16
to cusp-users
Steve,

I forget to tell you my cusp version in my previous post. It is v.0.5.1 which is the latest version. Please see my previous post for the errors.

Thanks.

Hongzhu

On Friday, July 22, 2016 at 10:18:57 AM UTC-6, Steve wrote:
Reply all
Reply to author
Forward
0 new messages