示例

/*******************************************************************************
    KML_EigenSolver example codes for solving eigenvalues problem AX = BXD 
*******************************************************************************/

#include "kes.h"

void MyMatVec4A(void *kesA, void *kesX, void *kesAX, int nRows, int nCols, int ld)
{
    // user-defined mat-vec routine for matrix A ......
}

void MyMatVec4B(void *kesB, void *kesX, void *kesBX, int nRows, int nCols, int ld)
{
    // user-defined mat-vec routine for matrix B ......
}

void MyMatVec4P(void *kesP, void *kesX, void *kesPX, int nRows, int nCols, int ld)
{
    // user-defined mat-vec routine for preconditioner P ......
}

int main (int argc, char **argv) {
    MPI_Init(&argc, &argv);
    int rank, size;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);

/*
    set parameters
*/
    int nRows = 1000, nCols = 100;
    MPI_Comm rowComm = MPI_COMM_WORLD;
    MPI_Comm colComm = MPI_COMM_SELF;

/*
    create matrix, vector, and solver
*/
    KesMatrix A, P, B;
    KesVector X, D;
    KesSolver solver;

    KesScalarType scalarType = KES_SCALAR_FP64;
    KesSolverAlgorithm solverAlgorithm = KES_LOBPCG;
    KesMatrixFormat matrixFormat = KES_SHELL;
    KesVectorFormat vectorFormat = KES_DENSE;

    KesMatrixCreate(&A, scalarType, matrixFormat);
    KesMatrixCreate(&B, scalarType, matrixFormat);
    KesMatrixCreate(&P, scalarType, matrixFormat);
    KesVectorCreate(&X, scalarType, vectorFormat);
    KesVectorCreate(&D, scalarType, vectorFormat);
    KesSolverCreate(&solver, scalarType, solverAlgorithm);

    KesMatrixSetDim(A, nRows * size, nRows * size);
    KesMatrixSetDim(B, nRows * size, nRows * size);

    KesVectorSetDim(X, nRows, nCols);
    KesVectorSetComm(X, KES_ROW, rowComm);
    KesVectorSetComm(X, KES_COL, colComm);

    KesVectorSetDim(D, nCols, 1);
    KesVectorSetComm(D, KES_ROW, rowComm);
    KesVectorSetComm(D, KES_COL, colComm);

    double *X_vals = new double[nRows * nCols];
    double *D_vals = new double[nCols];
    KesVectorDenseSetValues(X, X_vals, nRows);
    KesVectorDenseSetValues(D, D_vals, nCols);

    double *A_vals = new double[nRows * nRows];
    double *B_vals = new double[nRows * nRows];
    double *P_vals = new double[nRows];
    KesMatrixSetShellMatVec(A, A_vals, MyMatVec4A);
    KesMatrixSetShellMatVec(B, B_vals, MyMatVec4B);
    KesMatrixSetShellMatVec(P, P_vals, MyMatVec4P);

/*
    set solver-related parameters
*/
    int maxIter = 100000;
    int maxOrthoIter = 2;
    double threshold = 1e-6;
    KesProblemType problemType = KES_STANDARD;
    KesDenseEigenSolverType denseEigenSolverType = KES_SCALAPACK;
    KesOrthoProjectorType orthoProjectorType = KES_GRAM_SCHMIDT;
    KesOrthogonalizerType orthogonalizerType = KES_CHOLESKY_QR;
    KesConvergeDetectorType convergeDetectorType = KES_RESIDUALS;

    KesSolverSetThreshold(solver, threshold);
    KesSolverSetMaxIter(solver, maxIter);
    KesSolverSetMaxOrthoIter(solver, maxOrthoIter);
    KesSolverSetProblemType(solver, problemType);
    KesSetDenseEigenSolverType(solver, denseEigenSolverType);
    KesSetConvergeDetectorType(solver, convergeDetectorType);
    KesSetOrthogonalizerType(solver, orthogonalizerType);
    KesSetOrthoProjectorType(solver, orthoProjectorType);

/*
    run the solver
*/
    KesRunSolver(solver, A, B, P, X, D);

/*
    get solver information
*/
    int nIter;
    int nConv;
    KesGetSolverConv(solver, &nConv);
    KesGetSolverIter(solver, &nIter);
    if (rank == 0) {
        printf(" Number of iteration steps = %d\n", nIter);
        printf(" Number of converged eigenvalues = %d\n", nConv);
    }

/*
    compute residuals
*/
    KesVector R;
    KesVectorReplicate(X, &R, true, true);
    KesVector normR;
    KesVectorReplicate(D, &normR, true, true);
    KesEigResidual(A, B, X, D, R);
    KesVectorNorm(R, normR);

/*
    clean up
*/
    delete[] A_vals;
    delete[] B_vals;
    delete[] P_vals;
    delete[] X_vals;
    delete[] D_vals;

/*
    clean up internal object
*/
    KesMatrixDestroy(A);
    KesMatrixDestroy(B);
    KesMatrixDestroy(P);
    KesVectorDestroy(X);
    KesVectorDestroy(D);
    KesVectorDestroy(R);
    KesVectorDestroy(normR);
    KesSolverDestroy(solver);

    MPI_Finalize();
    return 0;
}