過程の計算処理を表示させたい。

Question

環境

• windows10
• visual studio 2017

/*********************************************
 * 連立方程式の解法 ( ガウスの消去法 )
 *********************************************/
#include <iostream>  // for cout
#include <stdio.h>   // for printf()

// 元の数定義
#define N 4  // 3

using namespace std;

/*
 * 計算クラス
 */
class Calc
{
    double a[N][N + 1];

    // 各種変数
    double d;     // ダミー
    int i, j, k;  // LOOP インデックス

    public:
        // 連立方程式を解く（ガウスの消去法）
        void calcGaussElimination();
};

/*
 * 連立方程式を解く（ガウスの消去法）
 */
void Calc::calcGaussElimination()
{
    // 係数
    static double a[N][N + 1] = {
        //{ 2.0, -3.0,  1.0,  5.0},
        //{ 1.0,  1.0, -1.0,  2.0},
        //{ 3.0,  5.0, -7.0,  0.0}
        { 1.0, -2.0,  3.0, -4.0,  5.0},
        {-2.0,  5.0,  8.0, -3.0,  9.0},
        { 5.0,  4.0,  7.0,  1.0, -1.0},
        { 9.0,  7.0,  3.0,  5.0,  4.0}
    };

    // 元の連立方程式をコンソール出力
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++)
            printf("%+fx%d ", a[i][j], j + 1);
        printf("= %+f\n", a[i][N]);
    }

    // 前進消去
    for (k = 0; k < N -1; k++) {
        for (i = k + 1; i < N; i++) {
            d = a[i][k] / a[k][k];
            for (j = k + 1; j <= N; j++)
                a[i][j] -= a[k][j] * d;
        }
    }

    // 後退代入
    for (i = N - 1; i >= 0; i--) {
        d = a[i][N];
        for (j = i + 1; j < N; j++)
            d -= a[i][j] * a[j][N];
        a[i][N] = d / a[i][i];
    }

    // 結果出力
    for (k = 0; k < N; k++)
        printf("x%d = %f\n", k + 1, a[k][N]);
}

/*
 * メイン処理
 */
int main()
{
    try
    {
        // 計算クラスインスタンス化
        Calc objCalc;

        // 連立方程式を解く（ガウスの消去法）
        objCalc.calcGaussElimination();
    }
    catch (...) {
        cout << "例外発生！" << endl;
        return -1;
    }

    // 正常終了
    return 0;
}

の以下の

 // 前進消去
    for (k = 0; k < N -1; k++) {
        for (i = k + 1; i < N; i++) {
            d = a[i][k] / a[k][k];
            for (j = k + 1; j <= N; j++)
                a[i][j] -= a[k][j] * d;
        }
    }

の部分において過程の計算処理が見たかったのでprintf関数を使ったのですが警告がでます。
変数の型を考慮しながら以下のように書きました。

// 前進消去
	for (k = 0; k < N - 1; k++) {
		for (i = k + 1; i < N; i++) {
			d = a[i][k] / a[k][k];
			for (j = k + 1; j <= N; j++)
				a[i][j] -= a[k][j] * d; //a[i][j] = a[i][j]-a[k][j] * d

			printf("a=%lf\nである", a);//付け加えた部分～～
			printf("%d=a[%d][%d] /a[%d][%d]\n ", a,i,k,k,k);
			printf("[k]=%d\n", k);
			printf("a[%d][%d]=a[%d][%d]-a[%d][%d]*%lf\n", i, j, i, j, k, j, d);//～～ここまで付け加えました。
		}
	}

以下が警告です。

1>------ ビルド開始: プロジェクト: サイトに載っていた連立方程式の計算, 構成: Release x64 ------
1>Source.cpp
1>c:\users\daito\source\repos\サイトに載っていた連立方程式の計算\サイトに載っていた連立方程式の計算\source.cpp(58): warning C4477: 'printf' : 書式文字列 '%lf' には、型 'double' の引数が必要ですが、可変個引数 1 は型 'double (*)[5]' です
1>c:\users\daito\source\repos\サイトに載っていた連立方程式の計算\サイトに載っていた連立方程式の計算\source.cpp(59): warning C4477: 'printf' : 書式文字列 '%d' には、型 'int' の引数が必要ですが、可変個引数 1 は型 'double (*)[5]' です
1>c:\users\daito\source\repos\サイトに載っていた連立方程式の計算\サイトに載っていた連立方程式の計算\source.cpp(59): warning C4313: 'printf': 書式文字列内の '%d' が引数 1 ('double (*)[5]' 型) と競合しています。
1>コード生成しています。
1>0 of 16 functions ( 0.0%) were compiled, the rest were copied from previous compilation.
1>  0 functions were new in current compilation
1>  0 functions had inline decision re-evaluated but remain unchanged
1>コード生成が終了しました。
1>サイトに載っていた連立方程式の計算.vcxproj -> c:\users\daito\source\repos\サイトに載っていた連立方程式の計算\x64\Release\サイトに載っていた連立方程式の計算.exe
1>プロジェクト "サイトに載っていた連立方程式の計算.vcxproj" のビルドが終了しました。
========== ビルド: 1 正常終了、0 失敗、0 更新不要、0 スキップ ==========

Answer 1

for (k = 0; k < N - 1; k++) {
    for (i = k + 1; i < N; i++) {
        d = a[i][k] / a[k][k];
        for (j = k + 1; j <= N; j++)
            a[i][j] -= a[k][j] * d; //a[i][j] = a[i][j]-a[k][j] * d

        printf("a=%lf\nである", a);//付け加えた部分～～
        //                      ^ (1)
        printf("%d=a[%d][%d] /a[%d][%d]\n ", a,i,k,k,k);
        //      ^^                           ^ (2)
        printf("[k]=%d\n", k);
        printf("a[%d][%d]=a[%d][%d]-a[%d][%d]*%lf\n", i, j, i, j, k, j, d);//～～ここまで付け加えました。
    }
}

こちらですが、上でコメントした 2 個所に問題があります。

(1) a は二次元配列なので printf() では表示できない。

前の回答者さんが指摘されている通りです。
ここ、本当は a を全部表示したいのではなく a[i][j] を表示したいのではありませんか。

(2) a は二次元配列なので printf() では表示できない。

ここ、本当は a ではなく d を表示したいのではありませんか。ちなみに d は double なのですが、 %d でよいですか。
クラス内の宣言で d がダミーとコメントされているのが気になりますが、 d は重要な役割を果たしているはず...。

まずはエラーや警告を読んでみるところから

前の回答者さんもおっしゃっていますが、警告は読まれましたか。たとえばこれ。

書式文字列 '%lf' には、型 'double' の引数が必要ですが、可変個引数 1 は型 'double (*)[5]' です

とても親切な警告文だと思います。これを見ると %lf に対して二次元配列を渡していることが一目で分かります。たとえ double (*)[5] が初めなんのことやら分からなかったとしても、 指摘されている行数を見る と %lf で表示されるべき場所に a を渡していることが分かりますよね。では a はなんだったかと考えると、二次元配列です。
そこで自分が本当は何をやりたかったのかもう一度考えてください。二次元配列の中身を全て表示したかったのなら、 for ループを書くことになるでしょう。特定の一つの要素を表示したくて添字を付け忘れているのなら、添字を追加してください。

エラーや警告は初心者のうちは「どうせ読んでも分からない」と考えがちなのかもしれません。確かに、 (特に C++ は) コンパイラが複雑怪奇なエラーを出すこともあります。でも、どれだけ分からないエラーが出ても、必ずそこには エラーや警告が起きた行数が書かれています から、一度は自分の書いたコードのその行を見直してください。エラーが 言っていることが分からなくても、間違いを直す方法は分かる 、ということは多々あります。

どうかエラーや警告を読むことを諦めないでくださいね。 どんなに難解でも、少なくとも指摘された行は見直す 、です。

ちなみに、それでもどうしても直し方が分からない場合は質問して大丈夫ですよ。そういうときのための teratail ですからね。自分がエラーや警告を読んで考えたことも質問に盛り込めば、回答者さんも「ちゃんと考えてから聞いているんだな」という印象になって、気持ち良く回答してくださると思いますよ。

(補足)

エラーや警告は、前のエラーにひきずられてどばーーーーーっと出ることがあります。大量のエラーが出ても一歩も引かないでくださいね。まずは一番最初のエラーを見てみましょう。たいてい最初のエラーは本物なので、そのエラーを直して再度コンパイルするとがばっとエラーが消えたりします。

あまり大規模なものだと毎回再コンパイルするわけにはいかないかもしれませんが、練習程度の小さい規模のものなら、一番最初のエラーをつぶしたらとりあえず再コンパイルしてみるとよいでしょう。

直すとエラーが増えるパターンというのもありますが、同じところに同じエラーが出ていないなら、それは確実に一歩進んでいます。諦めないで、また一番最初のエラーから読んでいきましょう。

応援していますので、少しずつやってみてくださいね。

Comments

[carnage0216]

親切に回答して頂き本当にありがとうございます。
じっくり読ませてエラー、警告を解決していきます。

Answer 2

おフザケついでに CUDA(cuSOLVER)で解いてみた。
この程度ではちっとも速くないけど。

// cusolver.cpp
/* how to compile : */
  set CUDA_DIR="C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v9.2"
  set CUDA_INC=%CUDA_DIR%/include
  set CUDA_LIB=%CUDA_DIR%/lib/x64
  cl -EHsc -I%CUDA_INC% -wd4819 cusolver.cpp cudart.lib cusolver.lib -link -libpath:%CUDA_LIB%
*/

#include <cuda_runtime.h>
#include <cusolverDn.h> // dense LAPACK

#include <cassert>
#include <iostream>
using namespace std;

template<typename T>
inline size_t bytesof(unsigned int s) { return s * sizeof(T); }

template<typename T>
T* allocate(unsigned int size) {
  T* result = nullptr;
  cudaError_t status = cudaMalloc(&result, bytesof<T>(size));
  assert( status == cudaSuccess );
  return result;
}

int main() {
  cusolverStatus_t status;

  // dense LAPACK
  cusolverDnHandle_t handle;
  status = cusolverDnCreate(&handle);
  assert( status == CUSOLVER_STATUS_SUCCESS );

  const int n = 4;

  float A[n][n] = {
    { 1.0, -2.0,  3.0, -4.0 },
    {-2.0,  5.0,  8.0, -3.0 },
    { 5.0,  4.0,  7.0,  1.0 },
    { 9.0,  7.0,  3.0,  5.0 }
  };
  float* dA = allocate<float>(n*n);
  cudaMemcpy(dA, A, bytesof<float>(n*n), cudaMemcpyHostToDevice);
  int lda = 4;

  // 必要なバッファ量を求め、確保する
  int worksize;
  status = cusolverDnSgetrf_bufferSize(
             handle,
             n,   // 行
             n,   // 列
             dA,  // A
             lda, // Aのヨコハバ
             &worksize);
  assert( status == CUSOLVER_STATUS_SUCCESS );
#ifdef _DEBUG
  cout << "worksize = " << worksize << endl;
#endif

  float* workspace = allocate<float>(worksize);

  // 計算結果に関する情報
  int* devInfo = allocate<int>(1);
  // ピボット
  int* pivot = allocate<int>(n);

  // LU分解 : dAに結果が求まる(それとpivot)
  status = cusolverDnSgetrf(
             handle,
             n,   // 行
             n,   // 列
             dA,  // A
             lda, // Aのヨコハバ
             workspace,
             pivot,
             devInfo);
#ifdef _DEBUG
  int info;
  cudaMemcpy(&info, devInfo, sizeof(int), cudaMemcpyDeviceToHost);
  cout << "info = " << info << endl;
#endif
  assert( status == CUSOLVER_STATUS_SUCCESS );

  float B[] = { 5.0, 9.0, -1.0, 4.0 };
  int nrhs = 1; // 問題数
  float* dB = allocate<float>(n*nrhs);
  cudaMemcpy(dB, B, bytesof<float>(n*nrhs), cudaMemcpyHostToDevice);
  int ldb = 4;

  // AX = B を解く (解XはBをoverrideする)
  status = cusolverDnSgetrs(
             handle,
             CUBLAS_OP_T,
             n,     // 行(=列)
             nrhs,  // 問題数
             dA,    // A
             lda,   // Aのヨコハバ
             pivot, // LU分解で得られたピボット
             dB,    // B
             ldb,   // Bのヨコハバ
             devInfo);
#ifdef _DEBUG
  cudaMemcpy(&info, devInfo, sizeof(int), cudaMemcpyDeviceToHost);
  cout << "info = " << info << endl;
#endif
  assert( status == CUSOLVER_STATUS_SUCCESS );

  // 結果を取得し、出力する
  cudaMemcpy(B, dB, bytesof<float>(n*nrhs), cudaMemcpyDeviceToHost);
  for ( int i = 0; i < n; ++i ) {
    cout << "x" << i+1 << " = " << B[i] << endl;
  }

  cudaFree(workspace);
  cudaFree(dA);
  cudaFree(dB);
  cudaFree(devInfo);
  cudaFree(pivot);

  cusolverDnDestroy(handle);
}

Comments

[carnage0216]

なんかここまでして頂いて本当にどうもありがとうございます。
ただ連立方程式の計算過程はパターン化しているのでGPUでのほうが早いと思ったのですが、もしかしたら計算以外の部分で時間を食ってしまいCPUと同等あるいは遅い結果になったのでしょうね。

[carnage0216]

正直に言いますと何をやっているのかはわかっているのですが、なんでこのように書けるのか頭が全く追いつきません。
ただただGPUをここまで使いこなせるのが羨ましい、凄いという感情だけが起きました。

Answer 3

解決したのでおマケ。classの体裁を整えた(+pivot操作を追加)。

/*********************************************
 * 連立方程式の解法 ( ガウスの消去法 )
 *********************************************/
#include <cstdio>   // printf()
#include <cmath>    // fabs()

#define PIVOT 1

/*
 * 計算クラス
 */
template<size_t N>
class Calc {
  double a[N][N + 1];

public:
  Calc(double c[N][N]);
  // ax = b を解く
  // 解xはb[N]に上書きする。
  bool calcGaussElimination(double b[N]);
};

template<size_t N>
Calc<N>::Calc(double c[N][N]) {
  for ( int i = 0; i < N; ++i ) {
    for ( int j = 0; j < N; ++j ) {
      a[i][j] = c[i][j];
    }
  }
}

/*
 * 連立方程式を解く（ガウスの消去法）
 */
template<size_t N>
bool Calc<N>::calcGaussElimination(double b[N]) {
  for ( int i = 0; i < N; ++i ) { a[i][N] = b[i]; }

  // 前進消去
  for ( int k = 0; k < N -1; k++ ) {
#if PIVOT
    double maxval = fabs(a[k][k]);
      int    maxidx = k;
      for ( int i = k+1; i < N; ++i ) {
        if ( maxval < fabs(a[i][k]) ) {
           maxval = fabs(a[i][k]);
           maxidx = i;
        }
      }
      if ( maxval == 0.0 ) return false;
      if ( maxidx != k ) {
        for ( int i = 0; i <= N; ++i ) {
          double t = a[maxidx][i];
          a[maxidx][i] = a[k][i];
          a[k][i] = t;
        }
      }
#endif
      for ( int i = k + 1; i < N; i++ ) {
        if ( a[k][k] == 0 ) {
          return false;
        }
        double d = a[i][k] / a[k][k];
        for (int j = k + 1; j <= N; j++ )
          a[i][j] -= a[k][j] * d;
      }
  }

  // 後退代入
  for ( int i = N - 1; i >= 0; i-- ) {
    if ( a[i][i] == 0 ) return false;
    double d = a[i][N];
    for (int j = i + 1; j < N; j++) {
      d -= a[i][j] * a[j][N];
    }
    a[i][N] = d / a[i][i];
  }

  // 結果出力
  for (int k = 0; k < N; k++)
    b[k] = a[k][N];

  return true;
}

/*
 * メイン処理
 */
int main() {

  const size_t N = 4;

  // 係数
  double a[N][N] = {
        { 1.0, -2.0,  3.0, -4.0 },
        {-2.0,  5.0,  8.0, -3.0 },
        { 5.0,  4.0,  7.0,  1.0 },
        { 9.0,  7.0,  3.0,  5.0 }
  };
  double b[N] = { 5.0, 9.0, -1.0, 4.0 };

  // 元の連立方程式をコンソール出力
  for (int i = 0; i < N; i++) {
    for (int j = 0; j < N; j++) {
      printf("%+fx%d ", a[i][j], j + 1);
    }
    printf("= %+f\n", b[i]);
  }
  printf("\n");

  // 計算クラスインスタンス化
  Calc<N> objCalc(a);

  // 連立方程式を解く（ガウスの消去法）
  if ( objCalc.calcGaussElimination(b) ) {
    for (int k = 0; k < N; k++) {
      printf("x%d = %f\n", k + 1, b[k]);
    }
  }

}

Answer 4

二次元配列のaに対して
printf("a=%lf\nである", a)
で出力されるわけないでしょ。
printfの何を勉強してるの？
てか、
printf("%+fx%d ", a[i][j], j + 1);
を
わざわざprintf("a=%lf\nである", a)に変更して二次元配列aがすべて表示されるとでも思ったとしか考えようがない。
修正前と後及びエラーメッセージとにらめっこして自身で対処しなされ。

Comments

[carnage0216]

＞＞修正前と後及びエラーメッセージとにらめっこして自身で対処しなされ。
そうさせて頂きます。
いつもいつもありがとうございます。