PETSc

    #include "petsc.h"

    /* 自定义帮助信息 */
    static char help[] = "Solves a tridiagonal linear system with KSP.\n\n";

    int
    main(int argc, char** args)
    {

        /* 初始化 */
        PetscErrorCode ierr;
        PetscMPIInt size;
        ierr = PetscInitialize(&argc, &args, (char*)0, help);
        if (ierr)
            return ierr;
        ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);
        CHKERRMPI(ierr);
        if (size != 1)
            SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

        /* 主程序 
         * ......
         */

        /* 退出 */
        ierr = PetscFinalize();
        return ierr;
    }

    Vec x;
    ierr = VecCreate(PETSC_COMM_WORLD, &x);
    CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject)x, "test");
    CHKERRQ(ierr);
    ierr = VecSetSizes(x, PETSC_DECIDE, H_len);
    CHKERRQ(ierr);
    ierr = VecSetFromOptions(x);
    CHKERRQ(ierr);
    Vec y;
    ierr = VecDuplicate(x, &y);

    int ix[1] = { 0 };
    double complex vx[1] = { 1 };
    ierr = VecSetValues(x, 1, ix, vx, INSERT_VALUES);
    CHKERRQ(ierr);
    ierr = VecAssemblyBegin(x);
    CHKERRQ(ierr);
    ierr = VecAssemblyEnd(x);
    CHKERRQ(ierr);

    struct timespec t0_mvm;
    timespec_get(&t0_mvm, TIME_UTC);

    ierr = MatMult(H, x, y);
    CHKERRQ(ierr);

    ierr = VecDestroy(&x);
    CHKERRQ(ierr);
    ierr = VecDestroy(&y);
    CHKERRQ(ierr);

    Mat A;
    ierr = MatCreate(PETSC_COMM_WORLD, &A);
    CHKERRQ(ierr);
    ierr = MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, A_len, A_len);
    CHKERRQ(ierr);
    ierr = MatSetFromOptions(A);
    CHKERRQ(ierr);
    ierr = MatSetUp(A);
    CHKERRQ(ierr);

    PetscScalar value[3] = {1,2,3};
    PetscInt col[3] = {0,1,2};
    PetscInt row[1] = {1};
    ierr = MatSetValues(A, 1, row, 2, col, value, INSERT_VALUES);
    CHKERRQ(ierr);
    ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    CHKERRQ(ierr);
    ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
    CHKERRQ(ierr);

    MatView(A,PETSC_VIEWER_STDOUT_WORLD);


    ierr = MatDestroy(&A);
    CHKERRQ(ierr);

    /* 定义每个矩阵都需要用到的一些数据 */
    typedef struct A_data {} A_data_t;

    /* 自定义的矩阵乘向量运算函数, y = A x */
    PetscErrorCode
    A_MatVec(Mat A, Vec x, Vec y)
    {

        /* 获取数据 */
        A_data_t* A_data;
        PetscErrorCode ierr;
        MatShellGetContext(A, &A_data);


        /* 获得本地的 x 向量数据 */
        PetscInt x_start, x_end;
        ierr = VecGetOwnershipRange(x, &x_start, &x_end);
        CAKERRQ(ierr);
        const PetscScalar* x_arr;
        ierr = VecGetArrayRead(x, &x_arr);
        CAKERRQ(ierr);

        /* y 向量要先置零 */
        ierr = VecSet(y, 0);
        CAKERRQ(ierr);

        /* 设置计算后 y 的值 */
        {
            PetscScalar new_value;
            ierr = VecSetValues(y, 1, &new_value, ADD_VALUES);
            CAKERRQ(ierr);
        }
        ierr = VecAssemblyBegin(y);
        CAKERRQ(ierr);
        ierr = VecAssemblyEnd(y);
        CAKERRQ(ierr);

        /* 结束计算要返还 x 的所有权 */
        ierr = VecRestoreArrayRead(x, &x_arr);
        CAKERRQ(ierr);
        return 0;
    }


    /* 创建一个临时向量，来获取 PETSc 自动分配的局域大小, 这里是方矩阵，长方形矩阵要分别获取两边的尺寸 */
    int A_local_size = 0;
    {
        Vec x;
        VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,A_len,&x);
        VecGetLocalSize(x,&A_local_size);
        VecDestroy(&x);
    }

    /* 创建无矩阵的矩阵 A，必须设置每个局域块 A 的大小 */
    Mat A;
    A_data_t A_data = {};
    ierr = MatCreateShell(PETSC_COMM_WORLD, A_local_size, A_local_size, A_len, A_len, &A_data, &A);
    CAKERRQ(ierr);

    /* 为 A 绑定矩阵乘向量运算 */
    ierr = MatShellSetOperation(Afree, MATOP_MULT, (void (*)(void))A_MatVec);
    CAKERRQ(ierr);