2012-03-29 12 views
0

私は正方形のnxn行列から上位のヘッセンベルグ行列に変換する関数(アッパー)を使って配列を送るプログラムを書いています。これは正常に実行されているが、私は問題がある別の関数を介して配列を送信する必要があります。配列の最初の要素と同じになるようにする必要があります。ただし、ダブル(mu)をゼロに設定しています。私はあなたのコードは動作しshoud配列のインデックスを取得する際の問題


#include <stdio.h> 
#include <math.h> 


int idx(int r, int c, int n) 
{ 
    return r * n + c; 
} 

void transpose(double *a, int n, double *at) 
{ 
    for (int i = 0; i < n; i++) { 
     for (int j = 0; j < n; j++) { 
      at[(idx (i, j, n))] = a[(idx (j, i, n))]; 
     } 
    } 
} 

void matrix_multiplication(double *a, double *b, int n, double *combo) 
{ 

    for (int i = 0; i < n; i++) { 
     for (int j = 0; j < n; j++) { 
      combo[(idx(i, j, n))] = 0; 
      for (int k = 0; k < n; k++) { 
       combo[(idx(i, j, n))] += a[(idx(i, k, n))] * b[(idx(k, j, n))]; 
      } 
     } 
    } 

} 

void identity(double *a, int n) 
{ 
    for (int i = 0; i < n; i++) { 
      for (int j = 0; j < n; j++) { 
       if (i == j) { 
        a[(idx(i, j, n))] = 1; 
       } 
       else { 
        a[(idx(i, j, n))] = 0; 
       } 
      } 
     } 
} 

void upperhes(double *a, int n, double *u, double *b) 
{ 

    //Sets b equal to a 
    for (int i = 0; i < (n*n); i++) { 
     b[i] = a[i]; 
    } 

    //Sets u equal to the identity matrix 
    identity(u, n); 

    int times = 0; 
    for (int i = 0; i < (n-2); i++) { 
     for (int j = n-1; j > (1 + times); j--) { 
      double c, s, r; 

      r = sqrt((b[(idx(j-1,i,n))] * b[(idx(j-1,i,n))]) + (b[(idx(j,i,n))] * b[(idx(j,i,n))])); 
      if (r < 0) { 
       r = (-1 * r); 
      } 
      if (r < pow(10,-50)) { 
       c = 1; 
       s = 0; 
      } 
      else { 
       c = b[(idx(j-1, i, n))]/r; 
       s = (-1 * b[(idx(j, i, n))])/r; 
      } 

      //store takes the original values of specific points in matrix b and stores them, so the values of b can be manipulated 
      double store[6]; 

      store[0] = b[(idx(j-1, i,n))]; 
      store[1] = b[(idx(j,i,n))]; 
      store[4] = b[(idx(i,j-1,n))]; 
      store[5] = b[(idx(i,j,n))]; 


      b[(idx (j,i,n))] = (s * store[0]) + (c * store[1]); 
      b[(idx(j-1,i,n))] = (c * store[0]) + (-s * store[1]); 

      for (int k= i+1; k<n; k++) { 

       store[2] = b[(idx(j-1,k,n))]; 
       store[3] = b[(idx(j,k,n))]; 

       b[(idx(j-1,k,n))] = (c * store[2]) - (s * store[3]); 
       b[(idx(j,k,n))] = (s * store[2]) + (c * store[3]); 
      } 


      b[(idx (i, j, n))] = (s * store[4]) + (c * store[5]); 
      b[(idx (i, j-1, n))] = (c * store[4]) - (s * store[5]); 

      for (int k = i + 1; k < n; k++) { 

       store[2] = b[(idx(k, j-1, n))]; 
       store[3] = b[(idx(k, j, n))]; 

       b[idx(k, j-1, n)] = (c * store[2]) - (s * store[3]); 
       b[idx(k, j, n)] = (s * store[2]) + (c * store[3]); 
      } 


      store[0] = u[(idx(j-1, 0, n))]; 
      store[1] = u[(idx(j, 0, n))]; 

      u[idx (j, 0, n)] = (s * store[0]) + (c * store[1]); 
      u[idx(j-1, 0, n)] = (c * store[0]) + (-s * store[1]); 

      for (int k= 1; k<n; k++) { 

       store[2] = u[idx(j-1, k, n)]; 
       store[3] = u[idx(j, k, n)]; 

       u[idx(j-1, k, n)] = (c * store[2]) - (s * store[3]); 
       u[idx(j, k, n)] = (s * store[2]) + (c * store[3]); 
      } 
     } 

     times++; 

    } 

    //Sets ut equal to the transpose of u 
    double ut[n*n]; 
    transpose(u, n, ut); 

    //Multiplies u and a together. 
    double ua[(n*n)]; 
    matrix_multiplication(u, a, n, ua); 

    double b_check[(n*n)]; 
    matrix_multiplication(ua, ut, n, b_check); 

    //Prints out the matrix! 
    for(int i=0; i<n; i++){ 
     for(int j =0; j<n; j++){ 
      printf("(%+.3f)\t", b[(idx (i,j,n))]); 
     } 
     printf("\n\n"); 
    } 
    printf("\n"); 

} 

void two_through_five(double *b, int n, double *eigenvalues, int total) 
{ 

    for (int times = 0; times < 100; times++){ 

     double mu = b[0]; 

     .......... 
    } 
} 

void qr_symmetric(double *a, int n, double *b) 
{ 

    //Creates a tridiagonal matrix by using a method that creates upper Hessenberg matrices on a symmetrical matrix a. 
    double u[n * n]; 
    upperhes(a, n, u, b); 

    //Creates an array to store the eigenvalues 
    double eigenvalues[n]; 

    two_through_five(b, n, eigenvalues, n); 
    for (int i = 0; i < n; i++) { 
     printf("%d\t", eigenvalues[n]); 
    } 
    printf("\n"); 

} 

int main(void) 
{ 
    int n; 
    n = 4; 
    double a[(n*n)]; 
    a[0] = 1; 
    a[1] = 2; 
    a[2] = 3; 
    a[3] = 4; 
    a[4] = 2; 
    a[5] = 6; 
    a[6] = 7; 
    a[7] = 8; 
    a[8] = 3; 
    a[9] = 7; 
    a[10] = 11; 
    a[11] = 15; 
    a[12] = 4; 
    a[13] = 8; 
    a[14] = 15; 
    a[15] = 16; 
    double b[n * n]; 
    qr_symmetric(a, n, b); 
    return 0; 
} 
+1

投稿したコードは、あなたがmuでゼロになっている理由の原因ではありません。あなたのupperhes()関数に何かがあるかもしれませんが、問題の短い、コンパイル可能な、完全な例を投稿する必要があります。 – tinman

+0

@tinman私の問題はもっと理にかなっているように私は残りのコードを掲載しました。 upperhessにはいくつかのバグがありますが、対称行列でうまくいくように見えますが、これは私の入力が確実なものです。 – GoodLuckGanesh

答えて

0

以下の私のコードが含まれていました。例えば、これは私の作品:

void two_through_five(double *b, int n, double *eigenvalues, int total) 
{ 
    int times; 
    for (times = 0; times < 100; times++) 
    { 
     double mu = b[0]; 
     printf("%f\n", mu); 
    } 
} 

int main() 
{ 

    int n = 18; 
    int i =0; 

    double * b = malloc(n*sizeof(double)); 

    for(i=0; i<n;i++) b[i] = 15.0; 

    double eigenvalues[n]; 

    two_through_five(b, n, eigenvalues, n); 

    return 0; 
} 

Outuptsを:

15.000000 
15.000000 
15.000000 
15.000000 
15.000000 
15.000000 
15.000000 
15.000000 
15.000000 
... 

だから私のquedtionは、あなたが[0] 0ではありませんBよろしいですか?

+0

私は残りのコードを投稿しました。このように実行すると、b [0]は1.000に等しくなります。 – GoodLuckGanesh

+0

two_through_five()を呼び出してからprintf( "%f \ n"、mu)を呼び出す直前にb [0]の値を出力しようとします。 inside two_through_five。多分それが価値を変えると思うかもしれません。なぜなら、私のコードがうまくいくとすれば、あなたのものもそうするべきだからです。私たちを掲載し続ける! – nay

+0

私はちょうど私の問題を認識しました。テストするために何かを印刷しているときは、printf( "%d"、b [0])またはprintf( "%d"、mu)を使っています。それは間違いなく動作します!ありがとう! – GoodLuckGanesh

関連する問題