C++ source code to do LU-decomposition by Crout method.

file.cpp

#include "file.hpp"

#include <iostream>
#include <sstream>
#include <string>
#include <vector>

bool File::get_text(std::vector<std::vector<double>>& data) {
  try {
    // ファイル OPEN
    std::ifstream ifs(f_data);
    if (!ifs.is_open()) return false;  // 読み込み失敗

    // ファイル READ
    std::string buf;                   // 1行分バッファ
    while (getline(ifs, buf)) {
      std::vector<double> rec;         // 1行分ベクタ
      std::istringstream iss(buf);     // 文字列ストリーム
      std::string field;               // 1列分文字列

      // 1行分文字列を1行分ベクタに追加
      double s;
      while (iss >> s)
        rec.push_back(s);

      // １行分ベクタを data ベクタに追加
      if (rec.size() != 0) data.push_back(rec);
    }
  } catch (...) {
      std::cerr << "EXCEPTION!" << std::endl;
      return false;
  }

  return true;  // 読み込み成功
}

file.hpp

#ifndef REGRESSION_LINE_FILE_HPP_
#define REGRESSION_LINE_FILE_HPP_

#include <fstream>
#include <string>
#include <vector>

class File {
  std::string f_data;

public:
  File(std::string f_data) : f_data(f_data) {}
  bool get_text(std::vector<std::vector<double>>&);
};

#endif

lu_decomposition.cpp

#include "lu_decomposition.hpp"

#include <iostream>
#include <sstream>
#include <vector>

/**
 * @brief      行列サイズ取得
 *
 * @return     行列サイズ(unsigned int)
 */
unsigned int LuDecomposition::get_size() {
  return data.size();
}

/**
 * @brief      LU 分解（クラウト法(Crout method)）
 *             （U の対角要素を 1 とする）
 *
 * @param[ref] m x m 行列（配列） mtx (double)
 * @return     真偽(bool)
 * @retval     true  成功
 * @retval     false 失敗
 */
bool LuDecomposition::lu_decompose_crout(
  std::vector<std::vector<double>>& mtx
) {
  std::size_t i;  // loop インデックス
  std::size_t j;  // loop インデックス
  std::size_t k;  // loop インデックス
  std::size_t m;  // 行列サイズ
  double s;       // 計算用(sum)
  double tmp;     // 一時退避用

  try {
    m = get_size();
    for (k = 0; k < m; k++) {
      for (i = k; i < m; i++) {
        s = 0.0;
        for (j = 0; j < k; j++)
          s += mtx[i][j] * mtx[j][k];
        mtx[i][k] -= s;
      }
      if (mtx[k][k] == 0) {
        std::cout << "[ERROR] Can't divide by 0!" << std::endl;
        return false;  // 計算失敗
      }
      tmp = 1.0 / mtx[k][k];
      for (j = k + 1; j < m; j++) {
        s = 0.0;
        for (i = 0; i < k; i++)
          s += mtx[k][i] * mtx[i][j];
        mtx[k][j] = (mtx[k][j] - s) * tmp;
      }
    }
  } catch (...) {
    return false;  // 計算失敗
  }

  return true;  // 計算成功
}

lu_decomposition.hpp

#ifndef REGRESSION_LINE_LU_DECOMPOSITION_HPP_
#define REGRESSION_LINE_LU_DECOMPOSITION_HPP_

#include <vector>

class LuDecomposition {
  std::vector<std::vector<double>> data;  // 元データ

public:
  LuDecomposition(std::vector<std::vector<double>>& data) : data(data) {}
                                          // コンストラクタ
  unsigned int get_size();                // 行列サイズ取得
  bool lu_decompose_crout(std::vector<std::vector<double>>&);
                                          // LU 分解（クラウト法）
};

#endif

lu_decomposition_crout.cpp

/***********************************************************
  LU 分解
  : クラウト法(Crout method)

    DATE          AUTHOR          VERSION
    2020.07.29    mk-mode.com     1.00 新規作成

  Copyright(C) 2020 mk-mode.com All Rights Reserved.
***********************************************************/
#include "lu_decomposition.hpp"
#include "file.hpp"

#include <cstdlib>   // for EXIT_XXXX
#include <iomanip>   // for setprecision
#include <iostream>
#include <string>
#include <vector>

int main(int argc, char* argv[]) {
  std::string f_data;                     // データファイル名
  std::vector<std::vector<double>> data;  // データ配列
  std::size_t i;                          // loop インデックス
  std::size_t j;                          // loop インデックス
  unsigned int m;                         // 行列サイズ

  try {
    // コマンドライン引数のチェック
    if (argc != 2) {
      std::cerr << "[ERROR] Number of arguments is wrong!\n"
                << "[USAGE] ./lu_decomposition_crout <file_name>"
                << std::endl;
      return EXIT_FAILURE;
    }

    // ファイル名取得
    f_data = argv[1];

    // データ取得
    File file(f_data);
    if (!file.get_text(data)) {
      std::cout << "[ERROR] Failed to read the file!" << std::endl;
      return EXIT_FAILURE;
    }

    // 計算用オプジェクトのインスタンス化
    LuDecomposition lu(data);

    // 行列サイズ
    m = lu.get_size();

    // データ一覧出力
    std::cout << "A = " << std::endl;
    std::cout << std::fixed << std::setprecision(4);
    for (i = 0; i < m; i++) {
      for (j = 0; j < m; j++)
        std::cout << std::setw(10) << std::right << data[i][j];
      std::cout << std::endl;
    }

    // 計算
    if (!lu.lu_decompose_crout(data)) {
      std::cout << "[ERROR] Failed to decompose!" << std::endl;
      return EXIT_FAILURE;
    }

    // 結果出力
    std::cout << "LU = " << std::endl;
    std::cout << std::fixed << std::setprecision(4);
    for (i = 0; i < m; i++) {
      for (j = 0; j < m; j++)
        std::cout << std::setw(10) << std::right << data[i][j];
      std::cout << std::endl;
    }
  } catch (...) {
      std::cerr << "EXCEPTION!" << std::endl;
      return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}