数组在通过 void 函数后没有改变

Question

我使用 Jacobi 算法用 C 编写了一个程序来计算矩阵的 SVD。我在程序中有 3 个函数：一个名为 jacobi_helper 的函数将执行 SVD，一个名为 jacobi 的函数在初始化几个数组后调用 jacobi_helper，以及调用 jacobi 的主函数。

我的问题在于，SVD 似乎在 jacobi_helper 中正确执行，但这些值似乎并不 return 到 jacobi，即使它们是在 void 函数中传递的.最后，我试图打印 SVD 的矩阵 U 和 V 以及奇异值向量，而我的困惑是由于奇异值向量能够很好地打印而矩阵却不能。

这是我的代码：

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

    double dot_product(double *a, double *b, int n){
        double sum = 0;
        for(int i = 0; i < n; i++){
            sum = sum + a[i]*b[i];
        }
        return sum;
    }

    //Right matrix must be transpose of the matrix you want to multiply
    void matrix_multiplication(left, right, output, n)
    int n;
    double left[n][n];
    double right[n][n];
    double output[n][n];
    {
        for(int i = 0; i < n; i++){
            for(int j = 0; j < n; j++){
                output[i][j] = dot_product(left[i], right[j], n);
            }
        }
    }

    void transpose(a, n, a_transpose)
    int n;
    double a[n][n];
    double a_transpose[n][n];
    {

        for(int i = 0; i < n; i++){
            for(int j = 0; j < n; j++){
                a_transpose[i][j] = a[j][i];
            }
        }
    }


    void jacobi_helper(a, n, givens, givens_transpose, s, s_matrix, u, temp_u, u_output, v, temp_v, v_output)
    int n; //Dimension of Matrix
    double a[n][n]; //Input Matrix
    double givens[n][n]; //Givens matrix to be formed in each step
    double givens_transpose[n][n]; //Transpose of Givens Matrix
    double s[n]; //Vector of singular values
    double s_matrix[n][n]; //Matrix containing the singular values.
    double u[n][n]; //U Matrix in SVD. (A = USV)
    double temp_u[n][n]; //Temporary Matrix
    double u_output[n][n];
    double v[n][n];
    double temp_v[n][n];
    double v_output[n][n];
    {
        u = a;


        //Set V to Identity Matrix
        for(int i = 0; i < n; i++){
            for(int j = 0; j < n; j++){
                if(i == j){
                    v[i][j] = 1;
                }
                else{
                    v[i][j] = 0;
                }
            }
        }

        double n_squared = 0;
        double value = 0;
        double epsilon = 0.1;

        for(int i = 0; i < n; i++){
            for(int j = 0; j < n; j++){
                n_squared = n_squared + (u[i][j]*u[i][j]);
            }
        }

        do{
            value = 0;
            for(int i = 0; i < (n-1); i++){
                for(int j = (i+1); j < n; j++){

                    double tmp1 = 0, tmp2 = 0, tmp4 = 0;
                    for(int k = 0; k < n; k++){
                        tmp1 = tmp1 + (u[k][i]*u[k][i]);
                        tmp2 = tmp2 + (u[k][i]*u[k][j]);
                        tmp4 = tmp4 + (u[k][j]*u[k][j]);
                    }

                    double tau = (tmp4 - tmp1)/(2*tmp2);

                    double t;
                    double t1;
                    double t2;

                    t1 = -tau + sqrt(1 + (tau*tau));
                    t2 = -tau - sqrt(1 + (tau*tau));

                    if(fabs(t1) < fabs(t2)){
                        t = t1;
                    }else{
                        t = t2;
                    }

                    double cos = 1/(sqrt(1+(t*t)));
                    double sin = cos*t;


                    //Make Givens Matrix
                    for(int k = 0; k < n; k++){
                        for(int l = 0; l < n; l++){
                            if(k == i && l == i){
                                givens[k][l] = cos;
                            }
                            else if (k == j && l == j){
                                givens[k][l] = cos;
                            }
                            else if ((k == i && l == j) && i < j){
                                givens[k][l] = -sin;
                            }
                            else if ((k == i && l == j) && i > j){
                                givens[k][l] = sin;
                            }
                            else if ((k == j && l == i) && i < j){
                                givens[k][l] = sin;
                            }
                            else if ((k == j && l == i) && i > j){
                                givens[k][l] = -sin;
                            }
                            else if(k == l){
                                givens[k][l] = 1;
                            }
                            else{
                                givens[k][l] = 0;
                            }
                        }
                    }


                    //Set U <- U*Givens
                    matrix_multiplication(u, givens, temp_u, n);
                    u = temp_u;

                    //Set V <- V*Givens
                    matrix_multiplication(v, givens, temp_v, n);
                    v = temp_v;

                    double temp = 0;
                    for(int k = 0; k < n; k++){
                        temp = temp + (u[k][i]*u[k][j]);
                    }
                    value = value + temp*temp;
                }
            }

        }while(sqrt(value) >= (epsilon*epsilon)*(n_squared));

        for(int i = 0; i < n; i++){
            double sing_value = 0;
            for(int k = 0; k < n; k++){
                sing_value = sing_value + (u[k][i]*u[k][i]);
            }
            s[i] = sqrt(sing_value);
            s_matrix[i][i] = 1/sqrt(sing_value);
        }

        matrix_multiplication(u, s_matrix, temp_u, n);
        u = temp_u;

        printf("%.16f  %.16f  %.16f  %.16f\n",u[0][0],u[0][1],u[1][0],u[1][1]);
        printf("%.16f  %.16f  %.16f  %.16f\n",v[0][0],v[0][1],v[1][0],v[1][1]);
        printf("%.16f  %.16f\n\n", s[0], s[1]);
    }


    void jacobi(double *a, int n, double *s, double *u, double *v){
        static double s_matrix[5000000];
        static double givens[5000000];
        static double givens_transpose[5000000];
        static double temp_u[5000000];
        static double temp_v[5000000];
        static double u_output[5000000];
        static double v_output[5000000];


        jacobi_helper(a, n, givens, givens_transpose, s, s_matrix, u, temp_u, u_output, v, temp_v, v_output);
    }

    int main(int argc, const char * argv[]) {
        static double a[4] = {2,2,-1,1};
        int n = 2;
        static double s[2];
        static double u[4];
        static double v[4];


        jacobi(a, n, s, u, v);
        printf("%.16f  %.16f  %.16f  %.16f\n",u[0],u[1],u[2],u[3]);
        printf("%.16f  %.16f  %.16f  %.16f\n",v[0],v[1],v[2],v[3]);
        printf("%.16f  %.16f\n", s[0], s[1]);
        return 0;
    }

所以问题是，当我在末尾有 printf 语句时，s[0] 和 s[1] 在两种情况下都是相同的，但是 u[0] 到 u[3] 和 v[0 ] 到 v[3] 是不一样的。

从 jacobi_helper 函数中的打印语句我得到（正确的值）：

    u = 1.0000000000000000  0.0000000000000000  0.0000000000000000  1.0000000000000000
    v = 0.7071067811865475  -0.7071067811865475  0.7071067811865475  0.7071067811865475
    s = 2.8284271247461898  1.4142135623730949

但是从主函数我得到：

    u = 0.0000000000000000  0.0000000000000000  0.0000000000000000  0.0000000000000000
    v = 1.0000000000000000  0.0000000000000000  0.0000000000000000  1.0000000000000000
    s = 2.8284271247461898  1.4142135623730949

有人可以帮忙吗？如果这有用，我使用的算法是来自该网站的算法 6：http://www.cs.utexas.edu/users/inderjit/public_papers/HLA_SVD.pdf

Answer 1

jacobi_helper 是用参数 double *u 调用的，它是指向存储结果的内存地址的指针，但 jacobi_helper 做的第一件事是 u = a 这意味着：忘记内存地址 u 并使用 a 代替。 您在指向参数的数组中设置一些值 v , 但后来你继续 v = temp_v.

数组变量只是一个指向内存地址的指针，为该数组变量赋值将用另一个地址替换该内存地址。如果你想改变数组变量指向的内存，你必须取消引用地址。

密码

u[0] = 123;
u[1] = 345;
u[2] = 678;

更改 u 指向的内存中保存的数据。另一方面，

u = a;

只会改变 u，所有对 u[0] 的后续访问实际上将访问与 a[0].

相同的内存

还有参数double *u和参数double u[n][n]的类型不匹配。

数组在通过 void 函数后没有改变

Array not changing after being passed through void function

c

arrays

matrix

linear-algebra

numerical-methods