示例

#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;
}