示例
C Interface：
#include <stdio.h>
#include <math.h>
#include <mpi.h>
#include <omp.h>
#include <assert.h>
#include "amg4c.h"
int test(int isonlyAMG)
{
int nthreads = 1;
int rtn = 0;
int gid;
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm, &gid);
int globalN = 8;
int localOffset = 0;
int localLength = globalN;
int offsets[9] = {0, 2, 4, 6, 7, 8, 10, 12, 14};
int indices[14] = {0, 7, 1, 6, 2, 5, 3, 4, 2, 5, 1, 6, 0, 7};
double values[14] = {1.0, 2.0, -2.0, 3.0, 3.0, 4.0, -4.0, 5.0, 4.0, -6.0, 3.0, 7.0, 2.0, 8.0};
double b_[8] = {3.0, 1.0, 7.0, -4.0, 5.0, -2.0, 10.0, 10.0};
double x_[8] = {0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5};
int isAvg = 0;
int64_t localNnz = offsets[localLength];
int *indices_ = NULL;
double *values_ = NULL;
indices_ = indices;
values_ = values;
KmlAmgMatrixH matrixA = NULL;
KmlAmgMatrixOptions MatrixAoptions;
KmlAmgMatrixStore MatrixAstore;
{
MatrixAoptions.fieldMask = KMLAMG_MATRIX_OPTIONS_VALUE_TYPE | KMLAMG_MATRIX_OPTIONS_INDEX_TYPE |
KMLAMG_MATRIX_OPTIONS_COMM | KMLAMG_MATRIX_OPTIONS_NUM_THREADS;
MatrixAoptions.comm = comm;
MatrixAoptions.indexType = KMLAMG_INDEX_INT32;
MatrixAoptions.valueType = KMLAMG_VALUE_FP64;
MatrixAoptions.numThreads = nthreads;
MatrixAstore.globalNRows = globalN;
MatrixAstore.globalNCols = globalN;
MatrixAstore.isAvgPart = isAvg;
MatrixAstore.localBegin = localOffset;
MatrixAstore.localNRows = localLength;
MatrixAstore.localNCols = localLength;
MatrixAstore.type = KMLAMG_MATRIX_STORE_CSR;
MatrixAstore.csr.offsets = offsets;
MatrixAstore.csr.indices = indices_;
MatrixAstore.csr.values = NULL;
if (KmlAmgMatrixCreate(&matrixA, &MatrixAstore, &MatrixAoptions) != KMLAMG_SUCCESS) {
printf("Create A error!\n");
};
}
if (KmlAmgMatrixSetValues(matrixA, values_) != KMLAMG_SUCCESS) {
printf("matrixA SetValues error!\n");
}
KmlAmgMatrixH VectorB = NULL;
KmlAmgMatrixOptions VectorBoptions;
KmlAmgMatrixStore VectorBstore;
{
VectorBoptions.fieldMask = KMLAMG_MATRIX_OPTIONS_VALUE_TYPE | KMLAMG_MATRIX_OPTIONS_INDEX_TYPE |
KMLAMG_MATRIX_OPTIONS_COMM | KMLAMG_MATRIX_OPTIONS_NUM_THREADS;
VectorBoptions.comm = comm;
VectorBoptions.indexType = KMLAMG_INDEX_INT32;
VectorBoptions.valueType = KMLAMG_VALUE_FP64;
VectorBoptions.numThreads = nthreads;
VectorBstore.globalNRows = globalN;
VectorBstore.globalNCols = 1;
VectorBstore.isAvgPart = isAvg;
VectorBstore.localBegin = localOffset;
VectorBstore.localNRows = localLength;
VectorBstore.localNCols = 1;
VectorBstore.type = KMLAMG_MATRIX_STORE_DENSE_ROW_MAJOR;
VectorBstore.dense.values = NULL;
if (KmlAmgMatrixCreate(&VectorB, &VectorBstore, &VectorBoptions) != KMLAMG_SUCCESS) {
printf("Create B error!\n");
};
}
if (KmlAmgMatrixSetValues(VectorB, b_) != KMLAMG_SUCCESS) {
printf("VectorB SetValues error!\n");
}
KmlAmgMatrixH VectorX = NULL;
KmlAmgMatrixOptions VectorXoptions;
KmlAmgMatrixStore VectorXstore;
{
VectorXoptions.fieldMask = KMLAMG_MATRIX_OPTIONS_VALUE_TYPE | KMLAMG_MATRIX_OPTIONS_INDEX_TYPE |
KMLAMG_MATRIX_OPTIONS_COMM | KMLAMG_MATRIX_OPTIONS_NUM_THREADS;
VectorXoptions.comm = comm;
VectorXoptions.indexType = KMLAMG_INDEX_INT32;
VectorXoptions.valueType = KMLAMG_VALUE_FP64;
VectorXoptions.numThreads = nthreads;
VectorXstore.globalNRows = globalN;
VectorXstore.globalNCols = 1;
VectorXstore.isAvgPart = isAvg;
VectorXstore.localBegin = localOffset;
VectorXstore.localNRows = localLength;
VectorXstore.localNCols = 1;
VectorXstore.type = KMLAMG_MATRIX_STORE_DENSE_ROW_MAJOR;
VectorXstore.dense.values = x_;
if (KmlAmgMatrixCreate(&VectorX, &VectorXstore, &VectorXoptions) != KMLAMG_SUCCESS) {
printf("Create X error!\n");
};
}
if (KmlAmgMatrixSetValues(VectorX, x_) != KMLAMG_SUCCESS) {
printf("VectorX SetValues error!\n");
}
KmlAmgPrecondH precond = NULL;
KmlAmgPrecondOptions PreOptions;
KmlAmgPrecondSetupOptions PreSetOptions;
{
PreOptions.fieldMask = KMLAMG_PRECOND_OPTIONS_VALUE_TYPE | KMLAMG_PRECOND_OPTIONS_INDEX_TYPE |
KMLAMG_PRECOND_OPTIONS_TYPE | KMLAMG_PRECOND_OPTIONS_NUM_THREADS;
PreOptions.valueType = KMLAMG_VALUE_FP64;
PreOptions.indexType = KMLAMG_INDEX_INT32;
PreOptions.precondType = KMLAMG_PRECOND_AMG;
PreOptions.numThreads = nthreads;
if (KmlAmgPrecondCreate(&precond, &PreOptions) != KMLAMG_SUCCESS) {
printf("Create Precond error!\n");
};
PreSetOptions.fieldMask =
KMLAMG_PRECOND_SETUP_OPTIONS_PARTITION | KMLAMG_PRECOND_SETUP_OPTIONS_STRONG_THRESHOLD |
KMLAMG_PRECOND_SETUP_OPTIONS_MAX_LEVELS | KMLAMG_PRECOND_SETUP_OPTIONS_COARSEN_TYPE |
KMLAMG_PRECOND_SETUP_OPTIONS_INTERP_TYPE | KMLAMG_PRECOND_SETUP_OPTIONS_DOWN_SMOOTH_TYPE |
KMLAMG_PRECOND_SETUP_OPTIONS_UP_SMOOTH_TYPE | KMLAMG_PRECOND_SETUP_OPTIONS_COARSEST_SMOOTH_TYPE;
PreSetOptions.strongThreshold = 0.25;
PreSetOptions.maxLevels = 1;
PreSetOptions.coarsenType = KMLAMG_PRECOND_AMG_COARSEN_PMIS;
PreSetOptions.interpType = KMLAMG_PRECOND_AMG_INTERP_DIRECT;
PreSetOptions.downSmoothType = KMLAMG_PRECOND_AMG_SMOOTH_L1_SOR_FORWARD;
PreSetOptions.upSmoothType = KMLAMG_PRECOND_AMG_SMOOTH_L1_SOR_BACKWARD;
PreSetOptions.coarsestSmoothType = KMLAMG_PRECOND_AMG_SMOOTH_L1_SOR_SYMMETRIC;
PreSetOptions.partitionAlgo = KMLAMG_PARTITION_KRDR;
if (KmlAmgPrecondSetSetupOptions(precond, &PreSetOptions) != KMLAMG_SUCCESS) {
printf("Set Precond setup options error!\n");
};
if (KmlAmgPrecondSetup(precond, matrixA) != KMLAMG_SUCCESS) {
printf("Setup Precond error!\n");
};
}
KmlAmgSolverH Ksp = NULL;
KmlAmgSolverOptions KspOptions;
KmlAmgSolverSetupOptions KspSetupOptions;
{
KspOptions.fieldMask = KMLAMG_SOLVER_OPTIONS_VALUE_TYPE | KMLAMG_SOLVER_OPTIONS_INDEX_TYPE |
KMLAMG_SOLVER_OPTIONS_TYPE | KMLAMG_SOLVER_OPTIONS_NUM_THREADS;
KspOptions.valueType = KMLAMG_VALUE_FP64;
KspOptions.indexType = KMLAMG_INDEX_INT32;
KspOptions.solverType = KMLAMG_SOLVER_CG;
KspOptions.numThreads = nthreads;
if (KmlAmgSolverCreate(&Ksp, &KspOptions) != KMLAMG_SUCCESS) {
printf("KmlAmgSolverCreate error!\n");
};
KspSetupOptions.fieldMask = KMLAMG_SOLVER_SETUP_OPTIONS_PARTITION;
KspSetupOptions.partitionAlgo = KMLAMG_PARTITION_KRDR;
}
if (isonlyAMG) {
KmlAmgPrecondApplyOptions preApplyOptions;
preApplyOptions.fieldMask = 0;
KmlAmgPrecondSetApplyOptions(precond, &preApplyOptions);
KmlAmgPrecondApply(precond, VectorB, VectorX);
KmlAmgPrecondInfo Precondinfo;
Precondinfo.fieldMask = 0;
KmlAmgPrecondQuery(precond, &Precondinfo);
double *xvalue = (double *)VectorXstore.dense.values;
printf("x is ");
for (int i = 0; i < globalN; i++) {
printf("%f ", xvalue[i]);
}
printf("\n");
return 0;
}
KmlAmgSolverSetPrec(Ksp, precond);
KmlAmgSolverApplyOptions KspApplyOptions;
{
KspApplyOptions.fieldMask = KMLAMG_SOLVER_APPLY_OPTIONS_MAX_ITERS | KMLAMG_SOLVER_APPLY_OPTIONS_TOLERANCE;
KspApplyOptions.maxIters = 1000;
KspApplyOptions.tolerance = 1e-9;
MPI_Barrier(comm);
for (int s = 0; s < 1; s++) {
KmlAmgMatrixSetValues(matrixA, values_);
KmlAmgMatrixSetValues(VectorB, b_);
MPI_Barrier(comm);
KmlAmgSolverSetSetupOptions(Ksp, &KspSetupOptions);
KmlAmgSolverSetup(Ksp, matrixA);
for (int ss = 0; ss < 1; ss++) {
MPI_Barrier(comm);
KmlAmgMatrixSetValues(VectorX, x_);
KmlAmgSolverSetApplyOptions(Ksp, &KspApplyOptions);
KmlAmgSolverApply(Ksp, VectorB, VectorX);
}
}
}
KmlAmgSolverInfo SolverInfo;
SolverInfo.fieldMask = 0;
KmlAmgSolverQuery(Ksp, &SolverInfo);
double *xvalue = (double *)VectorXstore.dense.values;
printf("x is ");
for (int i = 0; i < globalN; i++) {
printf("%f ", xvalue[i]);
}
printf("\n");
int *perm = (int *)malloc(sizeof(int)*globalN);
KmlAmgPrecondInfo precondInfo;
precondInfo.fieldMask = KMLAMG_PRECOND_INFO_LOCAL_BEGIN | KMLAMG_PRECOND_INFO_PERM;
precondInfo.perm = (void *)perm;
KmlAmgPrecondQuery(precond, &precondInfo);
int64_t permSum = 0;
for (int64_t i = 0; i < globalN; i++) {
permSum += perm[i];
}
assert(permSum == (globalN * (globalN - 1)) / 2);
KmlAmgSolverInfo solverInfo;
solverInfo.fieldMask = KMLAMG_SOLVER_INFO_LOCAL_BEGIN | KMLAMG_SOLVER_INFO_PERM | KMLAMG_SOLVER_INFO_NUM_ITERS;
solverInfo.perm = (void *)perm;
KmlAmgSolverQuery(Ksp, &solverInfo);
permSum = 0;
for (int64_t i = 0; i < globalN; i++) {
permSum += perm[i];
}
assert(permSum == (globalN * (globalN - 1)) / 2);
printf("nIterations = %ld\n", solverInfo.nIters);
KmlAmgMatrixDestroy(matrixA);
KmlAmgMatrixDestroy(VectorB);
KmlAmgMatrixDestroy(VectorX);
KmlAmgPrecondDestroy(precond);
KmlAmgSolverDestroy(Ksp);
return rtn;
}
int main(int argc, char **argv)
{
MPI_Init(&argc, &argv);
test(0);
MPI_Finalize();
return 0;
}
父主题： TYAMG