C Interface：#include "mg.h"
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include "mpi.h"
int main(int argc, char **argv)
{
MPI_Init(&argc, &argv);
MPI_Comm comm = MPI_COMM_WORLD;
int rankSize;
MPI_Comm_size(comm, &rankSize);
int pId;
MPI_Comm_rank(comm, &pId);
// Init proc size
int px = 1;
int py = 1;
int pz = 1;
// Init matrix size
int gx = 400;
int gy = 200;
int gz = 50;
int lx = gx / px;
int ly = gy / py;
int lz = gz;
int xbeg = (pId % px) * lx;
int xend = xbeg + lx;
int ybeg = (pId / px) * ly;
int yend = ybeg + ly;
int zbeg = 0;
int zend = lz;
int totLen = (lx + 2) * (ly + 2) * (lz + 2);
// Init matrix A
KmlMgMatrixH A = NULL;
float *aValue = (float *)malloc(19 * totLen * sizeof(float));
{
KmlMgMatrixStore store;
store.comm = comm;
store.globalx = gx;
store.globaly = gy;
store.globalz = gz;
store.procx = px;
store.procy = py;
store.procz = pz;
store.xbeg = xbeg;
store.xend = xend;
store.ybeg = ybeg;
store.yend = yend;
store.zbeg = zbeg;
store.zend = zend;
store.xhalo = 1;
store.yhalo = 1;
store.zhalo = 1;
store.crossPolar = 1;
store.type = KML_MG_MATRIX_STORE_SPARSE_HALO;
store.sparse.values = NULL;
store.sparse.type = KML_MG_SPARSE_AOS;
store.sparse.stencil = KML_MG_STENCIL_3D19;
KmlMgMatrixOptions options;
options.fieldMask = KML_MG_MATRIX_OPTIONS_INDEX_TYPE | KML_MG_MATRIX_OPTIONS_VALUE_TYPE;
options.indexType = KML_MG_INDEX_INT32;
options.valueType = KML_MG_VALUE_FP32;
if (A != NULL) {
KmlMgMatrixDestroy(A);
}
if (KmlMgMatrixCreate(&A, &store, &options) != KML_MG_OK) {
printf("Create A error!\n");
};
// Set A values, order is Y-X-Z
for (int i = 1; i < lx + 1; i++) {
for (int j = 1; j < ly + 1; j++) {
for (int k = 1; k < lz + 1; k++) {
for (int l = 0; l < 19; l++) {
aValue[j * (lx + 2) * (lz + 2) * 19 + i * (lz + 2) * 19 + k * 19 + l] = 0.5;
}
aValue[j * (lx + 2) * (lz + 2) * 19 + i * (lz + 2) * 19 + k * 19 + 9] *= 100;
}
}
}
if (KmlMgMatrixSetValues(A, (void *)aValue) != KML_MG_OK) {
printf("SetValues of A error!\n");
};
}
// Init B & X
KmlMgMatrixH B = NULL;
float *bValue = (float *)malloc(totLen * sizeof(float));
KmlMgMatrixH X = NULL;
float *xValue = (float *)malloc(totLen * sizeof(float));
{
KmlMgMatrixStore store;
store.comm = comm;
store.globalx = gx;
store.globaly = gy;
store.globalz = gz;
store.procx = px;
store.procy = py;
store.procz = pz;
store.xbeg = xbeg;
store.xend = xend;
store.ybeg = ybeg;
store.yend = yend;
store.zbeg = zbeg;
store.zend = zend;
store.xhalo = 1;
store.yhalo = 1;
store.zhalo = 1;
store.crossPolar = 1;
store.type = KML_MG_MATRIX_STORE_DENSE_COL_MAJOR_HALO;
store.dense.values = NULL;
store.dense.numDim = 1;
store.dense.is2D = 0;
store.dense.needCorner = 1;
KmlMgMatrixOptions options;
options.fieldMask = KML_MG_MATRIX_OPTIONS_INDEX_TYPE | KML_MG_MATRIX_OPTIONS_VALUE_TYPE;
options.indexType = KML_MG_INDEX_INT32;
options.valueType = KML_MG_VALUE_FP32;
// Create B
if (KmlMgMatrixCreate(&B, &store, &options) != KML_MG_OK) {
printf("Create B error!\n");
};
// Set B values
for (int i = 0; i < totLen; i++) {
bValue[i] = 1.1;
}
if (KmlMgMatrixSetValues(B, (void *)bValue) != KML_MG_OK) {
printf("SetValues of B error!\n");
};
// Create X
if (KmlMgMatrixCreate(&X, &store, &options) != KML_MG_OK) {
printf("Create X error!\n");
};
// Set X values
for (int i = 0; i < totLen; i++) {
xValue[i] = 0.0;
}
if (KmlMgMatrixSetValues(X, (void *)xValue) != KML_MG_OK) {
printf("SetValues of X error!\n");
};
}
// Init solver
KmlMgSolverH solver = NULL;
KmlMgPreconditionerOptions preOpt; // need malloc and free
{
KmlMgSolverOptions solOpt;
solOpt.fieldMask = KML_MG_SOLVER_OPTIONS_INDEX_TYPE | KML_MG_SOLVER_OPTIONS_KSP_VALUE_TYPE |
KML_MG_SOLVER_OPTIONS_PC_VALUE_TYPE | KML_MG_SOLVER_OPTIONS_PC_CALC_TYPE |
KML_MG_SOLVER_OPTIONS_KSP;
solOpt.indexType = KML_MG_INDEX_INT32;
solOpt.kspValueType = KML_MG_VALUE_FP32;
solOpt.pcValueType = KML_MG_VALUE_FP32;
solOpt.pcCalcType = KML_MG_VALUE_FP32;
solOpt.ksp = KML_MG_SOLVER_GCR;
if (solver != NULL) {
KmlMgSolverDestroy(solver);
}
if (KmlMgSolverCreate(&solver, comm, &solOpt) != KML_MG_OK) {
printf("solver Create error!\n");
};
preOpt.fieldMask = KML_MG_PRECONDITIONER_OPTIONS_INDEX_TYPE | KML_MG_PRECONDITIONER_OPTIONS_KSP_VALUE_TYPE |
KML_MG_PRECONDITIONER_OPTIONS_VALUE_TYPE | KML_MG_PRECONDITIONER_OPTIONS_CALC_TYPE |
KML_MG_PRECONDITIONER_OPTIONS_NUM_LEVEL | KML_MG_PRECONDITIONER_OPTIONS_SHIFT_LEVEL |
KML_MG_PRECONDITIONER_OPTIONS_CYCLE_TYPE | KML_MG_PRECONDITIONER_OPTIONS_SWEEP |
KML_MG_PRECONDITIONER_OPTIONS_XPSS | KML_MG_PRECONDITIONER_OPTIONS_SMOOTHER |
KML_MG_PRECONDITIONER_OPTIONS_WEIGHT | KML_MG_PRECONDITIONER_OPTIONS_COARSEN |
KML_MG_PRECONDITIONER_OPTIONS_RESTRICT | KML_MG_PRECONDITIONER_OPTIONS_PROLONG |
KML_MG_PRECONDITIONER_OPTIONS_ALPHA | KML_MG_PRECONDITIONER_OPTIONS_CONTROL;
preOpt.indexType = KML_MG_INDEX_INT32;
preOpt.kspValueType = KML_MG_VALUE_FP32;
preOpt.valueType = KML_MG_VALUE_FP32;
preOpt.calcType = KML_MG_VALUE_FP32;
preOpt.numLevels = 4;
preOpt.shiftLevelId = 1;
preOpt.cycle = KML_MG_CYCLE_V;
int *gridSweep = (int *)malloc(2 * sizeof(int));
gridSweep[0] = 1;
gridSweep[1] = 0;
preOpt.gridSweep = gridSweep;
int *xpss = (int *)malloc(4 * sizeof(int));
xpss[0] = 1;
xpss[1] = 1;
xpss[2] = 1;
xpss[3] = 1;
preOpt.xpss = xpss;
KmlMgSmootherType *smoother = (KmlMgSmootherType *)malloc(4 * sizeof(KmlMgSmootherType));
smoother[0] = KML_MG_SMOOTHER_PLANE_ILU;
smoother[1] = KML_MG_SMOOTHER_LINE_GS;
smoother[2] = KML_MG_SMOOTHER_LINE_GS;
smoother[3] = KML_MG_SMOOTHER_LINE_GS;
preOpt.smoother = smoother;
float *weight = (float *)malloc(4 * sizeof(float));
weight[0] = 1.0;
weight[1] = 1.0;
weight[2] = 1.0;
weight[3] = 1.0;
preOpt.weight = weight;
KmlMgCoarsenType *coarsen = (KmlMgCoarsenType *)malloc(4 * sizeof(KmlMgCoarsenType));
coarsen[0] = KML_MG_COARSEN_SEMI_XY;
coarsen[1] = KML_MG_COARSEN_SEMI_XY;
coarsen[2] = KML_MG_COARSEN_SEMI_XY;
coarsen[3] = KML_MG_COARSEN_SEMI_XY;
preOpt.coarsen = coarsen;
KmlMgRestInteType *rest = (KmlMgRestInteType *)malloc(4 * sizeof(KmlMgRestInteType));
rest[0] = KML_MG_REST_INTE_CELL_2D16;
rest[1] = KML_MG_REST_INTE_CELL_2D4;
rest[2] = KML_MG_REST_INTE_CELL_2D4;
rest[3] = KML_MG_REST_INTE_CELL_2D4;
preOpt.rest = rest;
KmlMgRestInteType *interp = (KmlMgRestInteType *)malloc(4 * sizeof(KmlMgRestInteType));
interp[0] = KML_MG_REST_INTE_CELL_2D16;
interp[1] = KML_MG_REST_INTE_CELL_2D16;
interp[2] = KML_MG_REST_INTE_CELL_2D16;
interp[3] = KML_MG_REST_INTE_CELL_2D16;
preOpt.interp = interp;
float *alpha = (float *)malloc(4 * sizeof(float));
alpha[0] = 1.0;
alpha[1] = 1.0;
alpha[2] = 1.0;
alpha[3] = 1.0;
preOpt.alpha = alpha;
int *control = (int *)malloc(4 * sizeof(int));
control[0] = 0;
control[1] = 1;
control[2] = 1;
control[3] = 1;
preOpt.control = control;
if (KmlMgSolverSetPreconditioner(solver, &preOpt) != KML_MG_OK) {
printf("solver SetPreconditioner error!\n");
}
printf("SetPreconditioner done!\n");
if (KmlMgSolverSetPreconditionerMatrix(solver, A) != KML_MG_OK) {
printf("solver SetPreconditionerMatrix error!\n");
}
printf("SetPreconditionerMatrix done!\n");
if (KmlMgSolverPrepare(solver, A) != KML_MG_OK) {
printf("solver Prepare error!\n");
}
printf("Prepare done!\n");
}
// Apply
{
KmlMgApplyOptions opt;
opt.fieldMask = KML_MG_APPLY_OPTIONS_USE_ZERO_GUESS | KML_MG_APPLY_OPTIONS_TOLERANCE_TYPE |
KML_MG_APPLY_OPTIONS_TOLERANCE | KML_MG_APPLY_OPTIONS_MAX_ITERATION |
KML_MG_APPLY_OPTIONS_RESTART_STEP;
opt.useZeroGuess = 1;
opt.toleranceType = KML_MG_TOLERANCE_ABS;
opt.tolerance = 1e-6;
opt.maxIteration = 200;
opt.restartStep = 10;
KmlMgStatus sta = KmlMgSolverApply(solver, B, X, &opt);
if (sta != KML_MG_OK) {
printf("solver Apply error!\n");
}
const char * str = KmlMgStatusStr(sta);
printf("%s\n", str);
printf("Apply done\n");
}
KmlMgMatrixInfo info;
info.fieldMask = KML_MG_MATRIX_INFO_STORE;
KmlMgMatrixQuery(X, &info);
// Clear
{
if (A != NULL) {
KmlMgMatrixDestroy(A);
}
if (aValue != NULL) {
free((void *)aValue);
}
if (B != NULL) {
KmlMgMatrixDestroy(B);
}
if (bValue != NULL) {
free((void *)bValue);
}
if (X != NULL) {
KmlMgMatrixDestroy(X);
}
if (xValue != NULL) {
free((void *)xValue);
}
if (preOpt.gridSweep != NULL) {
free((void *)preOpt.gridSweep);
}
if (preOpt.xpss != NULL) {
free((void *)preOpt.xpss);
}
if (preOpt.smoother != NULL) {
free((void *)preOpt.smoother);
}
if (preOpt.weight != NULL) {
free((void *)preOpt.weight);
}
if (preOpt.coarsen != NULL) {
free((void *)preOpt.coarsen);
}
if (preOpt.rest != NULL) {
free((void *)preOpt.rest);
}
if (preOpt.interp != NULL) {
free((void *)preOpt.interp);
}
if (preOpt.alpha != NULL) {
free((void *)preOpt.alpha);
}
if (preOpt.control != NULL) {
free((void *)preOpt.control);
}
}
MPI_Finalize();
return 0;
}