Portable, Extensible Toolkit for Scientific Computation/PETScの拡張とカスタマイズ

はじめに

PETScは、その拡張性と柔軟性が特徴の一つです。本章では、PETScを拡張して独自の機能を追加する方法や、既存の機能をカスタマイズする方法について解説します。具体的には、以下のトピックを取り上げます。

カスタムデータ構造の作成
カスタムソルバーの実装
既存ソルバーのカスタマイズ
外部ライブラリとの連携

カスタムデータ構造の作成

カスタムベクトルの作成

PETScでは、独自のベクトル型を作成して使用することができます。これにより、特定の問題に最適化されたデータ構造を利用できます。

例: カスタムベクトルの実装

#include "petscvec.h"

typedef struct {
    PetscScalar *data;
    PetscInt size;
} MyVector;

PetscErrorCode MyVecCreate(MPI_Comm comm, PetscInt n, Vec *v) {
    MyVector *myvec;
    PetscMalloc1(n, &myvec->data);
    myvec->size = n;
    VecCreateMPIWithArray(comm, 1, n, PETSC_DECIDE, myvec->data, v);
    return 0;
}

PetscErrorCode MyVecDestroy(Vec *v) {
    MyVector *myvec;
    VecGetArray(*v, (PetscScalar**)&myvec->data);
    PetscFree(myvec->data);
    VecDestroy(v);
    return 0;
}

int main(int argc, char **argv) {
    PetscInitialize(&argc, &argv, NULL, NULL);

    Vec v;
    MyVecCreate(PETSC_COMM_WORLD, 10, &v);

    // ベクトルの操作
    // ... (省略)

    MyVecDestroy(&v);
    PetscFinalize();
    return 0;
}

カスタム行列の作成

同様に、独自の行列型を作成することも可能です。

例: カスタム行列の実装

#include "petscmat.h"

typedef struct {
    PetscScalar **data;
    PetscInt rows, cols;
} MyMatrix;

PetscErrorCode MyMatCreate(MPI_Comm comm, PetscInt m, PetscInt n, Mat *A) {
    MyMatrix *mymat;
    PetscMalloc2(m, n, &mymat->data);
    mymat->rows = m;
    mymat->cols = n;
    MatCreateDense(comm, PETSC_DECIDE, PETSC_DECIDE, m, n, NULL, A);
    return 0;
}

PetscErrorCode MyMatDestroy(Mat *A) {
    MyMatrix *mymat;
    MatDenseGetArray(*A, (PetscScalar***)&mymat->data);
    PetscFree2(mymat->data);
    MatDestroy(A);
    return 0;
}

int main(int argc, char **argv) {
    PetscInitialize(&argc, &argv, NULL, NULL);

    Mat A;
    MyMatCreate(PETSC_COMM_WORLD, 10, 10, &A);

    // 行列の操作
    // ... (省略)

    MyMatDestroy(&A);
    PetscFinalize();
    return 0;
}

カスタムソルバーの実装

カスタム線形ソルバー

PETScでは、独自の線形ソルバーを実装して組み込むことができます。

例: カスタム線形ソルバーの実装

#include "petscksp.h"

PetscErrorCode MyLinearSolver(KSP ksp, Vec b, Vec x) {
    // カスタムソルバーの実装
    Vec r;
    VecDuplicate(b, &r);
    KSPGetOperators(ksp, &A, NULL);
    MatMult(A, x, r);
    VecAYPX(r, -1.0, b); // r = b - A*x
    VecCopy(r, x);
    VecDestroy(&r);
    return 0;
}

int main(int argc, char **argv) {
    PetscInitialize(&argc, &argv, NULL, NULL);

    KSP ksp;
    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, A, A);
    KSPSetType(ksp, KSPPREONLY); // 前処理のみ使用
    KSPSetComputeOperators(ksp, MyLinearSolver, NULL);
    KSPSetFromOptions(ksp);
    KSPSolve(ksp, b, x);

    PetscFinalize();
    return 0;
}

カスタム非線形ソルバー

同様に、非線形ソルバーもカスタマイズできます。

例: カスタム非線形ソルバーの実装

#include "petscsnes.h"

PetscErrorCode MyNonlinearSolver(SNES snes, Vec x, Vec f, void *ctx) {
    // カスタム非線形ソルバーの実装
    Vec r;
    VecDuplicate(f, &r);
    SNESComputeFunction(snes, x, r);
    VecAYPX(r, -1.0, f); // r = f - F(x)
    VecCopy(r, x);
    VecDestroy(&r);
    return 0;
}

int main(int argc, char **argv) {
    PetscInitialize(&argc, &argv, NULL, NULL);

    SNES snes;
    SNESCreate(PETSC_COMM_WORLD, &snes);
    SNESSetFunction(snes, r, FormFunction, NULL);
    SNESSetJacobian(snes, J, J, FormJacobian, NULL);
    SNESSetComputeFunction(snes, MyNonlinearSolver, NULL);
    SNESSetFromOptions(snes);
    SNESSolve(snes, NULL, x);

    PetscFinalize();
    return 0;
}

既存ソルバーのカスタマイズ

前処理のカスタマイズ

PETScでは、既存のソルバーに独自の前処理を追加することができます。

例: カスタム前処理の実装

#include "petscpc.h"

PetscErrorCode MyPreconditioner(PC pc, Vec x, Vec y) {
    // カスタム前処理の実装
    VecCopy(x, y);
    return 0;
}

int main(int argc, char **argv) {
    PetscInitialize(&argc, &argv, NULL, NULL);

    PC pc;
    PCCreate(PETSC_COMM_WORLD, &pc);
    PCSetType(pc, PCSHELL);
    PCShellSetApply(pc, MyPreconditioner);
    PCSetFromOptions(pc);

    // ソルバーの設定
    // ... (省略)

    PetscFinalize();
    return 0;
}

外部ライブラリとの連携

外部ライブラリの利用

PETScは、他の数値計算ライブラリ（例: SLEPc, Hypre, SuperLU）と連携することができます。

例: SLEPcとの連携

#include "petsc.h"
#include "slepceps.h"

int main(int argc, char **argv) {
    PetscInitialize(&argc, &argv, NULL, NULL);
    EPS eps;

    EPSCreate(PETSC_COMM_WORLD, &eps);
    EPSSetOperators(eps, A, NULL);
    EPSSetFromOptions(eps);
    EPSSolve(eps);

    // 固有値の取得
    // ... (省略)

    EPSDestroy(&eps);
    PetscFinalize();
    return 0;
}

練習問題

カスタムベクトル: 独自のベクトル型を作成し、簡単なベクトル演算を行うプログラムを作成してください。
カスタムソルバー: 独自の線形ソルバーを実装し、簡単な線形方程式を解くプログラムを作成してください。
前処理のカスタマイズ: 既存のソルバーに独自の前処理を追加し、性能を評価してください。
外部ライブラリの利用: SLEPcを使用して、行列の固有値問題を解くプログラムを作成してください。

まとめ

本章では、PETScを拡張して独自の機能を追加する方法や、既存の機能をカスタマイズする方法について学びました。カスタムデータ構造の作成、カスタムソルバーの実装、既存ソルバーのカスタマイズ、外部ライブラリとの連携を通じて、PETScの柔軟性と拡張性を活用する方法を理解しました。