QR分解C程序

#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include "./QRdecomposition.h"


#define dim		3

// 初始化矩阵  
void initializeMatrix(double* matrix, int rows, int cols) 
{
	int i = 0; 
	for(i = 0; i < rows * cols; i++)
	{
		matrix[i] = i + 1;  // 用1到rows*cols的值初始化矩阵  
	}
}


// 打印矩阵  
void printMatrix(double* matrix, int rows, int cols)
{
	int i = 0;
	int j = 0; 
	for(i = 0; i < rows; i++)
	{
		for (j = 0; j < cols; j++)
		{
			printf("%.3f ", matrix[i * cols + j]);
		}
		printf("\n");
	}
}



int main()
{
	int res = 0;
	double A[dim][dim] = {0};
	double A1[dim][dim] = {0};

	double Q[dim][dim] = {0};
	double R[dim][dim] = {0};

	//初始化矩阵A 
	initializeMatrix(&A[0][0], dim, dim);
	
	//打印矩阵A
	printf("Martix A:\n");
	printMatrix(&A[0][0], dim, dim);
	printf("\n");

	//QR分解函数
	res = qrDecomposition(dim, &A[0][0], &Q[0][0], &R[0][0]);
	if(res != 0)
	{
		printf("error:input dim > MAX_SIZE!\n");
		return -1;	
	}

	//矩阵与矩阵相乘
	cw_multiply(dim, dim, dim, &Q[0][0], &R[0][0], &A1[0][0]);

	//打印矩阵Q 
	printf("Martix Q:\n");
	printMatrix(&Q[0][0], dim, dim);
	printf("\n");
	
	//打印矩阵R	
	printf("Martix R:\n");
	printMatrix(&R[0][0], dim, dim);
	printf("\n");	
	
	//打印结果A1 
	printf("Martix A1:\n");
	printMatrix(&A1[0][0], dim, dim);
	
	return 0;
}

#include <string.h>
#include <math.h>
#include "QRdecomposition.h"
/************************************************************************************************************************
*  版权所有（C）2024，重庆长安工业（集团）有限责任公司。
*  
*  文件名称 ：      QRdecomposition.c
*  内容摘要 ：      QR分解
*  其他说明 ：      无
*  当前版本 ：      V1.00
*  作    者 ：      yf、yhz
*  完成日期 ：      20241024
*  修改记录 ：      
*  修改日期 ：      
*  版 本 号 ：
*  修 改 人 ：      
*  修改内容 ：
*************************************************************************************************************************/

/*
 * description:
 * 生成单位矩阵
 * input：
 * int n            矩阵大小
 * output：
 * double *matrix   单位矩阵首地址
 * return：
 * 无
*/
void cw_identityMatrix(int n, double *matrix)
{
    int i = 0;
    int j = 0;
    for(i = 0; i < n; i++)
    {
        for(j = 0; j < n; j++)
        {
            matrix[i*n + j] = (i == j) ? 1.0 : 0.0;
        }
    }
}

/*
 * description:
 * 两矩阵相乘
 * input：
 * int n            矩阵1行数
 * int m            矩阵1列数
 * int p            矩阵2列数
 * double *matrix1  矩阵1首地址
 * double *matrix2  矩阵2首地址
 * output：
 * double *result   结果矩阵首地址
 * return：
 * 无
 */
void cw_multiply(int n, int m, int p, double *matrix1, double *matrix2, double *result)
{
    int i = 0;
    int j = 0;
    int k = 0;
    for(i = 0; i < n; i++)
    {
        for(j = 0; j < p; j++)
        {
            result[i*p + j] = 0;
            for(k = 0; k < m; k++)
            {
                result[i*p + j] += matrix1[i*m + k] * matrix2[k*p + j];
            }
        }
    }
}

/*
 * description:
 * 矢量和标量相乘
 * input：
 * double *a 矢量首地址
 * double  b 标量值
 * int     m 矢量长度
 * output：
 * double *c 结果矢量首地址
 * return：
 * int		0：正常   其他：异常
*/
int cw_vecMulsca(double *a, double b, int m, double *c)
{
    int i = 0;

    if(a == NULL)
    {
        return -1;
    }
    if(c == NULL)
    {
        return -2;
    }
    if(m < 0)
    {
        return -3;
    }

    for(i = 0; i < m; i++)
    {
        c[i] = a[i] * b;
    }

    return 0;
}


/*
 * description:
 * QR分解算法
 * input：
 * int     n 矩阵维数
 * double *A 矩阵A首地址
  * 
 * output：
 * double *Q 结果正交矩阵Q首地址
 * double *R 结果上三角矩阵R首地址
 * return：
 * int		0：正常   其他：异常
*/
int qrDecomposition(int n, double *A, double *Q, double *R)
{
	int i = 0;
	int j = 0;
	int k = 0;
	double norm = 0.0;
	double w[MAX_SIZE] = {0.0};
	double x[MAX_SIZE] = {0.0};
	double E[MAX_SIZE*MAX_SIZE] = {0.0};
	double F[MAX_SIZE*MAX_SIZE] = {0.0};
	double P[MAX_SIZE*MAX_SIZE] = {0.0};
	double P1[MAX_SIZE*MAX_SIZE] = {0.0};
	double Y[MAX_SIZE*MAX_SIZE] = {0.0};

	if(n > MAX_SIZE)
	{
		return -1;	
	}

	cw_identityMatrix(n, E);             //单位矩阵E
	cw_identityMatrix(n, P1);            //单位矩阵P1

	for(k = 0; k < (n-1); k++)
	{
		norm = 0;
		for(i = k; i < n; i++)
		{
			norm += A[i*n + k]*A[i*n + k];
		}
		norm = sqrt(norm);             //求向量的范数

		if(A[k*n + k] > 0)
		{
			norm = -norm;
		}

		R[k*n + k] = -norm;

		for(i = 0; i < n; i++)
		{
			if(i < k)
			{
				w[i] = 0;
				R[i*n + k] = A[i*n + k];
			}
			else if(i == k)
			{
				w[i] = A[k*n + k] + norm;
			}
			else
			{
				w[i] = A[i*n + k];
			}
		}

		norm = 0;
		for(i = 0; i < n; i++)
		{
			norm += w[i]*w[i];
		}

		norm = sqrt(norm);
		if(fabs(norm) > 1e-6)
		{
			for(i = 0; i < n; i++)
			{
				w[i] = w[i]/norm;             //向量归一化
			}
		}

		cw_vecMulsca(w, 2.0, n, x);          
		cw_multiply(n, 1, n, x, w, F);       

		for(i = 0; i < n; i++)
		{
			for(j = 0; j < n; j++)
			{
				P[i*n + j] = E[i*n + j] - F[i*n + j];
			}
		}

		cw_multiply(n, n, n, P, A, Y);
		memcpy(&A[0], &Y[0], n*n*sizeof(double));

		cw_multiply(n, n, n, P1, P, Y);
		memcpy(&P1[0], &Y[0], n*n*sizeof(double));
		memcpy(&Q[0], &P1[0], n*n*sizeof(double));

		for(i = 0; i < n; i++)
		{
			R[i*n + (n-1)] = A[i*n + (n-1)];
		}
	}
	
	return 0;
}