私は正方形の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;
}
投稿したコードは、あなたがmuでゼロになっている理由の原因ではありません。あなたのupperhes()関数に何かがあるかもしれませんが、問題の短い、コンパイル可能な、完全な例を投稿する必要があります。 – tinman
@tinman私の問題はもっと理にかなっているように私は残りのコードを掲載しました。 upperhessにはいくつかのバグがありますが、対称行列でうまくいくように見えますが、これは私の入力が確実なものです。 – GoodLuckGanesh