/example_sparse_multiple_rhs.c Secret
Created
January 11, 2018 13:04
Star
You must be signed in to star a gist
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// This is a simple standalone example. See README.txt | |
#include <stdlib.h> | |
#include <stdio.h> | |
#include <string.h> | |
#include <math.h> | |
#include "magma_v2.h" | |
#include "magmasparse.h" | |
// ------------------------------------------------------------ | |
// This is an example how magma can be integrated into another software. | |
int main( int argc, char** argv ) | |
{ | |
// The software does e.g. discretization of a PDE, | |
// ends up with a sparse linear system in CSR format and a RHS. | |
// Let's assume this system is a diagonal system of size m. | |
int i, m=3, n=3; | |
int *col, *row; | |
double *val, *rhs, *sol; | |
magma_imalloc_cpu(&row, (m+1)); | |
magma_imalloc_cpu(&col, m); | |
magma_dmalloc_cpu(&val, m); | |
magma_dmalloc_cpu(&rhs, m*n); | |
magma_dmalloc_cpu(&sol, m*n); | |
for (i = 0; i < m; ++i) { | |
col[i] = i; | |
row[i] = i; | |
val[i] = 55.0; | |
for (int j = 0; j < n; ++j) | |
{ | |
rhs[i + j * m] = j + 3; | |
sol[i + j * m] = j; | |
} | |
} | |
row[m] = m; | |
// Initialize MAGMA and create some LA structures. | |
magma_init(); | |
magma_dopts opts; | |
magma_queue_t queue; | |
magma_queue_create( 0, &queue ); | |
magma_d_matrix A={Magma_CSR}, dA={Magma_CSR}; | |
magma_d_matrix b={Magma_CSR}, db={Magma_CSR}; | |
magma_d_matrix x={Magma_CSR}, dx={Magma_CSR}; | |
// Pass the system to MAGMA. | |
magma_dcsrset( m, m, row, col, val, &A, queue ); | |
magma_dvset( m, n, rhs, &b, queue ); | |
magma_dvset( m, n, sol, &x, queue ); | |
// Choose a solver, preconditioner, etc. - see documentation for options. | |
opts.solver_par.solver = Magma_CG; | |
opts.solver_par.restart = 8; | |
opts.solver_par.maxiter = 1000; | |
opts.solver_par.rtol = 1e-10; | |
opts.solver_par.maxiter = 1000; | |
opts.precond_par.solver = Magma_ILU; | |
opts.precond_par.levels = 0; | |
opts.precond_par.trisolver = Magma_CUSOLVE; | |
// Initialize the solver. | |
magma_dsolverinfo_init( &opts.solver_par, &opts.precond_par, queue ); | |
// Copy the system to the device (optional, only necessary if using the GPU) | |
magma_dmtransfer( A, &dA, Magma_CPU, Magma_DEV, queue ); | |
magma_dmtransfer( b, &db, Magma_CPU, Magma_DEV, queue ); | |
magma_dmtransfer( x, &dx, Magma_CPU, Magma_DEV, queue ); | |
// Generate the preconditioner. | |
magma_d_precondsetup( dA, db, &opts.solver_par, &opts.precond_par, queue ); | |
// In case we only wanted to generate a preconditioner, we are done. | |
// The preconditioner in the opts.precond_par structure - in this case an ILU. | |
// The lower ILU(0) factor is in opts.precond_par.L and | |
// the upper ILU(0) factor is in opts.precond_par.U (in this case on the device). | |
// If we want to solve the problem, run: | |
int code = magma_d_solver( dA, db, &dx, &opts, queue ); | |
if (code != MAGMA_SUCCESS) | |
{ | |
printf("Solver returned error code %d\n", code); | |
} | |
// Then copy the solution back to the host... | |
magma_dmfree( &x, queue ); | |
magma_dmtransfer( dx, &x, Magma_DEV, Magma_CPU, queue ); | |
// and back to the application code | |
magma_dvget( x, &m, &n, &sol, queue ); | |
// Free the allocated memory... | |
magma_dmfree( &dx, queue ); | |
magma_dmfree( &db, queue ); | |
magma_dmfree( &dA, queue ); | |
// and finalize MAGMA. | |
magma_queue_destroy( queue ); | |
magma_finalize(); | |
// From here on, the application code may continue with the solution in sol... | |
for (i = 0; i < m; ++i) { | |
for (int j = 0; j < n; ++j) | |
{ | |
printf("%.4f ", sol[i + j * m]); | |
} | |
printf("\n"); | |
} | |
return 0; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment