KmlScaissCgSolve?I

Solve a system of linear algebra equations. Each process reads and stores the global solution result.

Interface Definition

C interfaces:

int KmlScaissCgSolveSI(KmlScasolverTask **handle, int nb, float *x, int ldx, const float *b, int ldb);

int KmlScaissCgSolveDI(KmlScasolverTask **handle, int nb, double *x, int ldx, const double *b, int ldb);

Parameters

Parameter	Type	Description	Input/Output
handle	KmlScasolverTask **	Solver handle, which is a variable passed from a previous step.	Input/Output
nb	int	Number of RHS vectors. The value of nb can only be 1.	Input
x	float * for KmlScaissCgSolveSI double * for KmlScaissCgSolveDI	Array that reads and stores the global solution result. The number of elements n is the size of the global matrix.	Input/Output
ldx	int	Size of the leading dimension of x. This is a reserved parameter.	Input
b	const float * for KmlScaissCgSolveSI const double * for KmlScaissCgSolveDI	Local RHS array corresponding to the current process. The number of elements is the number of rows in the matrix corresponding to the current process, that is, stripeWidth.	Input
ldb	int	Size of the leading dimension of b. This is a reserved parameter.	Input

Return Value

Return Value	Type	Description
KMLSS_NO_ERROR	int	The execution is successful.
KMLSS_BAD_NB	int	The value of nb is not 1.
KMLSS_NO_MEMORY	int	The memory is insufficient.
KMLSS_NULL_ARGUMENT	int	A null parameter exists in {handle, x, b}.

Dependencies

#include "kml_scaiss.h"

Example

C interface:

MPI_Init(NULL, NULL);
int size, rank;
MPI_Comm_size(MPI_COMM_WORLD, &size);// Obtain the total number of processes.
MPI_Comm_rank(MPI_COMM_WORLD, &rank);// Obtain the current process rank.

int mat_size = 8;
int n = mat_size / size;
int n_beg = n * rank;
if (n * size != mat_size && rank == (size - 1)) {
    n = mat_size - n * rank;
}
int ja[26] = { 0, 3, 4, 1, 2, 3, 5, 1, 2, 7, 0, 1, 3, 6, 0, 4, 5, 1, 4, 5, 7, 3, 6, 2, 5, 7 };
double a[26] = { 1.0, 1.0, 2.0, 9.0, 2.0, 1.0, -3.0, 2.0, 3.0, 2.0, 1.0, 1.0, 9.0, -5.0, 2.0, 6.0, 1.0, -3.0, 1.0, 4.0, 1.0, -5.0, 7.0, 2.0, 1.0, 2.0 };
int ia[9] = { 0, 3, 7, 10, 14, 17, 21, 23, 26 };
int a_beg = ia[n_beg];
for (int i = n_beg; i < (n_beg + n + 1); i++) {
    ia[i] -= a_beg;
}

/* Internal KML_SCAISS structure */
KmlScasolverTask *handle;

/* KML_SCAISS control parameters */
int error;            /* Output error handle */

/* Create data structures */
const double *a_holder = &a[a_beg];
const int *ja_holder = &ja[a_beg];
const int *ia_holder = &ia[n_beg];
error = KmlScaissCgInitStripesDI(&handle, mat_size, 1, &n, &n_beg, &a_holder, &ja_holder, &ia_holder, MPI_COMM_WORLD);

double eps = 1e-4;
error = KmlScaissCgSetDID(&handle, KMLSS_THRESHOLD, &eps, 1); 
if (error != 0) {
    printf("\nERROR in KmlScaissCgSetDID with KMLSS_THRESHOLD: %d", error);
    return 1;
}
int max_iters = 2000;
error = KmlScaissCgSetDII(&handle, KMLSS_MAX_ITERATION_COUNT, &max_iters, 1);
if (error != 0) {
    printf("\nERROR in KmlScaissCgSetDII with KMLSS_MAX_ITERATION_COUNT: %d", error);
    return 1;
}
int precond = KMLSS_BJACOBI;
error = KmlScaissCgSetDII(&handle, KMLSS_PRECONDITIONER_TYPE, &precond, 1);
if (error != 0) {
    printf("\nERROR in KmlScaissCgSetDII with KMLSS_PRECONDITIONER_TYPE: %d", error);
    return 1;
}
int num_blocks = 1;
error = KmlScaissCgPcSetDII(&handle, KMLSS_NUM_BLOCKS, &num_blocks, 1);
if (error != 0) {
	printf("ERROR in KmlScaissCgPcSetDII: %d\n", error);
	return 1;
}
int subprecond = KMLSS_ICC;
error = KmlScaissCgPcSetDII(&handle, KMLSS_BLOCK_METHOD, &subprecond, 1);
if (error != 0) {
    printf("ERROR in KmlScaissCgPcSetDII: %d\n", error);
    return 1;
}
error = KmlScaissCgAnalyzeDI(&handle);
if (error != 0) {
    printf("ERROR in KmlScaissCgAnalyzeDI: %d\n", error);
    return 1;
}
error = KmlScaissCgFactorizeDI(&handle);
if (error != 0) {
    printf("ERROR in KmlScaissCgFactorizeDI: %d\n", error);
    return 1;
}

int nrhs = 1;
int ldx = mat_size, ldb = n;
double b[8] = { 4.0, 9.0, 7.0, 6.0, 9.0, 3.0, 2.0, 5.0 };
double x[8] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
error = KmlScaissCgSolveDxDI(&handle, nrhs, &x[0], ldx, &b[n_beg], ldb);
if (error != 0) {
    printf("ERROR in KmlScaissCgSolveDI: %d\n", error);
    return 1;
}

Parent topic: SCAISS-CG