scanfで値を読み込む所に、代わりに計算結果を代入したいのですが。。

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら　　※コードを貼り付ける場合は [code] と [/code] で囲って下さい。詳しくはこちら
閉鎖
yu

scanfで値を読み込む所に、代わりに計算結果を代入したいのですが。。

#1

投稿記事 by yu » 15年前

下のJacob法で行列の固有値、固有ベクトルを求める際に、求めたい行列Aの成分をscanfで読み取っているのですが、そこの部分を計算で求めた答えを直接代入する方法はないでしょうか？
例えば、

コード: 

/*Sの計算*/
 for(r=1; r<=n; r++)
	for(s=1; s<=n; s++)
	   for(i=1; i<=3; i++)
		  for(j=1; j<=3; j++)
    S[r][s]+=N[r][i]*N[s][j]*d[r][i]*d[s][j]*pow(PI/(a[r][i]+a[s][j]),1.5)
                                        *exp(-a[r][i]*a[s][j]*sqr(R[r][s])/(a[r][i]+a[s][j]));
のような式で求めたSという2×2行列の成分を代入したいのですが、、どなたか助けてください。


コード: 

//matrix.h(簡易行列操作ライブラリ)
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAXROW      20  /* 行数 */
#define MAXCOL      20  /* 列数 */
#define NUMITEMS    10  /* 画面に表示する際の1行の要素 */
#define NEARLY_ZERO    1.E-16 /* ゼロ */
#define ZERO_TOLERANCE 1.E-6  /* ゼロ */
#define vect_squared(n,x) vect_prod((n),(x),(x))
#define L2_norm(n,x) sqrt(vect_squared((n),(x)))

コード: 

/* Jacobi法 */
#define MAXROW      20  /* 行数 */
#define MAXCOL      20  /* 列数 */
#define NUMITEMS    10  /* 画面に表示する際の1行の要素 */
#define NEARLY_ZERO    1.E-16 /* ゼロ */
#define ZERO_TOLERANCE 1.E-6  /* ゼロ */
#define vect_squared(n,x) vect_prod((n),(x),(x))
#define L2_norm(n,x) sqrt(vect_squared((n),(x)))
typedef struct
{
	double R;
	double I;
} complex_t;
complex_t determinant(int n, complex_t A[][MAXCOL]);
double vect_prod(int n, double *x, double *y);
void mat_vect(int m, int n, double A[][MAXCOL], double *x, double *b);
void vect_mat(int m, int n, double *x, double A[][MAXCOL], double *b);
void mat_prod(int m, int n, int l,
			  double A[][MAXCOL],  double B[][MAXCOL], double C[][MAXCOL]);
void vect_copy(double *b, double *a, int n);
void mat_copy(int m, int n, double A[][MAXCOL], double B[][MAXCOL]);
void mat_read(int m, int n, double A[][MAXCOL]);
void vect_print(char *s, int n, double *x);
void mat_print(char *s, int m, int n, double A[][MAXCOL]);
double cleanup(int n, double A[][MAXCOL]);
#include "matrix.h"
/* ベクトルの内積 */
double vect_prod(int n, double *x, double *y)
{
	int i;
	double prod = 0.0;
	for(i = 0; i < n; i++)
		prod += x[i] * y[i];
	return prod;
}

/* 行列とベクトルの積 */
void mat_vect(int m, int n, double A[][MAXCOL], double *x, double *b)
{
	int i, j;
	for(i = 0; i < m; i++)
	{
		b[i] = 0.0;
		for(j = 0; j < n; j++)
			b[i] += A[i][j] * x[j];
	}
}
/* ベクトルと行列の積 */
void vect_mat(int m, int n, double *x, double A[][MAXCOL], double *b)
{
	int i, j;
	for(j = 0; j < n; j++)
	{
		b[j] = 0.0;
		for(i = 0; i < m; i++)
			b[j] += x[i] * A[i][j];
	}
}

/* 行列と行列の積 */
void mat_prod(int m, int n, int l,
			  double A[][MAXCOL], double B[][MAXCOL], double C[][MAXCOL])
{
	int i;
	for(i = 0; i < m; i++)
		vect_mat(n, n, &A[i][0], B, &C[i][0]) ;
}

/* ベクトルのコピー */
void vect_copy(double *b, double *a, int n)
{
	int i;
	for(i = 0; i < n; i++)
		b[i] = a[i];
}

/* 行列のコピー */
void mat_copy(int m, int n, double A[][MAXCOL], double B[][MAXCOL])
{
	int i;
	for(i = 0; i < m; i++)
		vect_copy(&B[i][0], &A[i][0], n);
}

/* 行列の読み込み */
void mat_read(int m, int n, double A[][MAXCOL])
{
	int i, j;
	for(i = 0; i < m; i++)
		for(j = 0; j < n; j++)
			scanf("%lf", &A[i][j]);
}

/* ベクトルの表示 */
void vect_print(char *s, int n, double *x)
{
	int j;
	if(s) printf("(%s)\n", s);
	for(j = 0; j < n; j++)
		printf("%11.4E%s", x[j], (j+1)% NUMITEMS ? " " : "\n");
	if(j % NUMITEMS)
		printf("\n");
}

/* 行列の表示 */
void mat_print(char *s, int m, int n, double A[][MAXCOL])
{
	int i;
	printf("(%s)\n", s);
	for(i = 0; i < m; i++)
		vect_print(NULL, n, &A[i][0]);
}

/* 微小要素の消去 */
double cleanup(int n, double A[][MAXCOL])
{
	int i, j;
	double  max_norm = 0.0;
	for(i = 0; i < n; i++)
		for(j = 0; j < n; j++)
			if(fabs(A[i][j]) > max_norm) max_norm = fabs(A[i][j]);
	if(max_norm > NEARLY_ZERO)
	{
		for(i = 0; i < n; i++)
			for(j = 0; j < n; j++)
				if(fabs(A[i][j]) / max_norm < ZERO_TOLERANCE)
					A[i][j] = 0.0;
	}
	return max_norm;
} 

//行列の固有値(古典的Jacob法)

#define MAXITR 200 /* 繰り返し回数の上限　*/

/* 回転行列 */
typedef struct
{
 int p, q; /* ゼロにする非対角要素の位置　*/
 double c; /* cos θ */
 double s; /* sin θ */
 double t; /* tan θ */
 double t2; /* s/(c+1) */
}JROTATE_t;

/* Aの更新 */
void update_A(int n, double A[][MAXCOL], JROTATE_t *J)
{
 int i;
 double aip, aiq;
 for(i = 0; i < n; i++)
 {
  aip = A[i][J->p];
  aiq = A[i][J->q];
  if(i != J->p && i != J->q)
  {
   A[i][J->p] = aip - J->s * (J->t2 * aip + aiq);
   A[i][J->q] = aiq + J->s * (aip - J->t2 * aiq);
   A[J->p][i] = A[i][J->p];
   A[J->q][i] = A[i][J->q];
  }
  else if(i == J->p)
  {
   A[i][i] = aip - J->t * aiq;
   A[i][J->q] = 0.0;
  }
  else
  {
   A[i][i] = aiq + J->t * aip;
   A[i][J->p] = 0.0;
  }
 }
}

/* Vの更新 */
void update_V(int n, double V[][MAXCOL], JROTATE_t *J)
{
 int i;
 double vip, viq;
 for(i = 0; i < n; i++)
 {
  vip = V[i][J->p];
  viq = V[i][J->q];
  V[i][J->p] = vip - J->s * (J->t2 * vip + viq);
  V[i][J->q] = viq + J->s * (vip - J->t2 * viq);
 }
}

/* 次にゼロにする非対角要素の選択(古典的ヤコビ法)*/
void select_nondiag(int n, double A[][MAXCOL], JROTATE_t *J)
{
 if(J->p == n - 2) /* 最終行*/
 {
  J->p = 0;
  J->q = 1;
 }
 else if(J->q == n - 1) /* 最終列*/
 {
  J->p++;
  J->q = J->p + 1;
 }
 else
    J->q++;
}

/* 回転行列を求める*/
void calc_Jmat(double A[][MAXCOL], JROTATE_t *J)
{
 if(fabs(A[J->p][J->q]) < NEARLY_ZERO)
 {
  J->c = 1.0;
  J->s = 0.0;
  J->t2 = 0.0;
 }
 else
 {
  double phi;
  phi = 0.5 * (A[J->q][J->q] - A[J->p][J->p]) / A[J->p][J->q];
  if(fabs(phi) < 1.E10)
      J->t = 1.0 / (fabs(phi) + sqrt(phi*phi + 1.0));
  else
     J->t = 0.5 / fabs(phi);
  J->t = phi < 0.0 ?  -J->t : J->t;
  J->c = 1.0 / sqrt(1.0 + J->t * J->t);
  J->s = J->t * J->c;
  J->t2 = J->s / (1.0 + J->c);
 }
}

/* 非対角要素の２乗和*/
double nondiag_norm(int n, double A[][MAXCOL])
{
 int i, j;
 double norm = 0.0;
 for(i = 0; i < n; i++)
 {
  for(j = 0; j < n; j++)
  {
   if(j == i)
      continue;
   norm += A[i][j] * A[i][j];
  }
 } 
 return sqrt(norm);
}

/*
* ヤコビ法
* 実対称行列の固有値、固有ベクトルの計算
*　　(入力)
*　　　　n : 行列サイズ
*　　　　A[][] : nxn実対称行列
*　　(出力)
*　　　　lamda[] : Aの固有値
*　　　　A[][] : 固有ベクトルを並べたもの
*　　(戻り値)
*　　　　0 : 正常に終了
*　　　　1 : 繰り返し回数の上限に達した
*/

int Jacobi(int n, double A[][MAXCOL], double lambda[])
{
 int i, j, iteration = 0, ret_code;
 double V[MAXROW][MAXCOL];
 JROTATE_t J;
	J.p = 0;
	J.q = 0;
 for(i = 0; i < n; i++)
 {
	 for(j = 0; j < n; j++)
	 V[i][j] = 0.0;
     V[i][i] = 1.0;
 }
 mat_print("A", n, n, A);
 mat_print("V", n, n, V);
 while(1)
 {
  if(nondiag_norm(n, A) < ZERO_TOLERANCE)
  {
   ret_code = 0; break; /* 正常終了 */
  }
  if(iteration++ >= MAXITR)
  {
   ret_code = 1; break; /* 繰り返し回数の上限 */
  }
  printf("\n----- Iteration %d -----\n", iteration);
  select_nondiag(n, A, &J);
  calc_Jmat(A, &J);
  printf("(p,q) = (%d, %d) c = %f s = %f\n", J.p, J.q, J.c, J.s);
  update_A(n, A, &J);
  update_V(n, V, &J);
 cleanup(n, A); cleanup(n, V);
  mat_print("A", n, n, A);
  mat_print("V", n, n, V);
  printf("norm of nondiagonals = %f\n", nondiag_norm(n, A));
 }
 for(i = 0; i < n; i++)
  lambda[i] = A[i][i];
 mat_copy(n, n, V, A);
 return ret_code;
}

/* メイン関数 */
int main(void)
{
 int n, ret_status;
 double A[MAXROW][MAXCOL], Aorg[MAXROW][MAXCOL];
 double lambda[MAXCOL];
 printf("n×n正方行列のnを指定してください。\n");
 printf("n:");  scanf("%d", &n);
 mat_read(n, n, A);
 mat_print("Input", n, n, A);
 mat_copy(n, n, A, Aorg);
 ret_status = Jacobi(n, A, lambda);
 vect_print("eigen values", n, lambda);
 printf("status = %d\n", ret_status);
 return 0;
}
ページトップ
yu

Re: scanfで値を読み込む所に、代わりに計算結果を代入したいのですが。。

#2

投稿記事 by yu » 15年前

メイン関数内の、mat_readの部分どうにかすれば、できるようになるのでしょうか？
自分はまだCを始めたばかりでどう修正すれば良いか良く解りません。よろしくお願いします。
ページトップ
non
記事: 1097
登録日時: 15年前

Re: scanfで値を読み込む所に、代わりに計算結果を代入したいのですが。。

#3

投稿記事 by non » 15年前

プログラムは見てもわかりませんので見ませんが、
最初のｓの計算部分が、本体のJacobi法と独立しているのでしょうか？
質問の意味がよく理解できていないのですが、最初の配列Sが求まるのなら、これをファイルに出力し、Jacobiのプログラムのほうでは、そのファイルをリダイレクトで読み込めばいいのではないでしょうか。
non
ページトップ
yu

Re: scanfで値を読み込む所に、代わりに計算結果を代入したいのですが。。

#4

投稿記事 by yu » 15年前

non さんが書きました:プログラムは見てもわかりませんので見ませんが、
最初のｓの計算部分が、本体のJacobi法と独立しているのでしょうか？
質問の意味がよく理解できていないのですが、最初の配列Sが求まるのなら、これをファイルに出力し、Jacobiのプログラムのほうでは、そのファイルをリダイレクトで読み込めばいいのではないでしょうか。
fprintf,fscanfなどを使えば良いのでしょうか？試して見たのですが、計算結果をDATA(ファイル名)に出力はできたのですが、その値をJacobiのプログラムの方へ読み込ませることができませんでした。まだ初めて3ヶ月くらいで理解が足らない部分が多々あるので、もう少し詳しく上の方法を教えていただけないでしょうか?お手数おかけます。
ページトップ
閉鎖

“C言語何でも質問掲示板” へ戻る