SVD分解（奇异值分解）的C语言实现

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

double determinant(double matrix[20][20],int order);
double laplace_expansion(double matrix[20][20], int r, int c,int order);


int max_row_order(int n, double A[n][n])
{
    double m = 0;
    int mro = 0;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (j == i)
            {
                continue;
            }
            
            if (fabs(A[i][j]) > m)
            {
                m = fabs(A[i][j]);
                mro = i;
            }  
        }                 
    }
    int p = mro;
    return p;
}

int max_column_order(int n, double A[n][n])
{
    double m = 0;
    int mco = 0;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (j == i)
            {
                continue;
            }
            
            if (fabs(A[i][j]) > m)
            {
                m = fabs(A[i][j]);
                mco = j;
            }  
        }                 
    }
    int q = mco;
    return q;
}

double t_slove(int n, double A[n][n], int p, int q)
{
    double s = (A[q][q] - A[p][p])/(2*A[p][q]);
    double t;
    if (s == 0)
    {
        t = 1;
    }
    else
    {
        double x1 = -s + sqrt(s*s+1);
        double x2 = -s - sqrt(s*s+1);
        if (fabs(x1) > fabs(x2))
        {
            t = x2;
        }
        else
        {
            t = x1;
        }
    }
    return t;
}

double calculate_nondiagonal(int n, double A[n][n])
{
    double sum = 0;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (j == i)
            {
                continue;
            }
            sum = sum + A[i][j] * A[i][j];
        }                 
    }
    return sum;
}

void print_matrix(int row, int col, double matrix[row][col])
{
    for(int i = 0;i < row;i ++)
    {
        for(int j = 0;j < col;j ++)
        {
            printf("%.12lf\t",matrix[i][j]);
        }    
        printf("\n");
    }
    printf("----------------------------------------------");
    printf("\n");
}

double calculate_mode(int row, double x[row][1])
{
    double mode2 = 0;
    for (int i = 0; i < row; i++)
    {
        mode2 = mode2 + x[i][0]*x[i][0];
    }
    double mode = sqrt(mode2);
    return mode;
}

double * GaussEliminationWithPartialPivoting1(int n, double A[n-1][n-1], double b[n-1])
{
    double *result;
    int size = n-1;
    result = (double *)malloc(size*sizeof(double));     
    for (int i = 0; i < n-1; i++)
    {
        result[i] = 0;
    }

    for (int i = 0; i < n-1; i++)
    {
        int k = i;
        for (int j = i+1; j < n-1; j++)
        {
            if (fabs(A[k][i]) < fabs(A[j][i]))
            {
                k = j;
            }
        }
        for (int m = 0; m < n-1; m++)
        {
            double c;
            c = A[i][m];
            A[i][m] = A[k][m];
            A[k][m] = c;
        }
        double d;
        d = b[i];
        b[i] = b[k];
        b[k] = d;

        for (int j = i+1; j < n-1; j++)
        {
            double v = -A[j][i]/A[i][i];
            for (int k = 0; k < n-1; k++)
            {
                A[j][k] = v*A[i][k] + A[j][k];
            }
            b[j] = v*b[i] + b[j];
        }  
    }