QR分解求矩阵逆--c工程实现

一、逆上三角矩阵算法
问题
我有代码来计算下三角矩阵的逆。如何通过稍微改变它来从下面的代码 计算上三角矩阵的逆？

function L = L_inv(A)
            [n,n] = size(A);
            L = zeros(n); 
            for i=1:n
                L(i,i) = 1/A(i,i);
                for j=i+1:n
                    L(j,i)=-A(j, i:j-1)*L(i:j-1,i)/A(j,j);
                end
            end
        end

笔者根据上述图片写的代码：

在这里插入代码片
/**
 * Function:Solution for Matrix Inverse 
 * Date:2022-07-16 18:00:00
 * Last Modified:Juchunyu
 */
#include <stdio.h>
#include <math.h>

`int main(){

  double m_in_[3][3]={
    1,2,3,
    0,2,3,
    0,0,5,
  };
  int length = 3;
  double m_in[3][3]={0};
  //tranpose
  for(int i=0;i<length;i++){
    for(int j=0;j<length;j++){
       m_in[j][i] = m_in_[i][j];
    }
  }
  
  double m_inverse_[3][3]={0};
  for(int j=0;j<length;j++){
      //求解对角线
      m_inverse_[j][j] = 1/m_in[j][j];
      for(int i=j+1;i<length;i++){
             double temp = 0;
             for(int k=j;k<=i-1;k++){
               temp+=m_in[i][k]*m_inverse_[k][j];
             }
             m_inverse_[i][j] = -temp/m_in[i][i];
      }
  }
 double m_inverse[3][3]={0};
//tranpose
  for(int i=0;i<length;i++){
    for(int j=0;j<length;j++){
       m_inverse[j][i] = m_inverse_[i][j];
    }
  }

  printf("shangsanjiao分解：\n");
  //printf
 for(int i=0;i<length;++i)
  {
      for(int j=0;j<length;++j)
      {
          printf("%2.4f    ",m_inverse[i][j]);
      //  
      }
      printf("\n");
  }

}`

二、QR分解

笔者根据上述的QR分解的数值计算公式写了如下的代码，可以求出Q和R矩阵。

/** * Function:Solution for Matrix Inverse * Date:2022-07-16 22:00:00 * Last Modified:Juchunyu */
#include <stdio.h>
#include <math.h>

#define SIZE 3


int main(){
    
    //double init[3][3]={1,1,1,1,1,0,1,0,0};
    double init[3][3]={
    1,2,3,4,5,6,7,8,9};
    double Q[3][3]={
    0};
    double R[3][3]={
    0};
    double v[3]={
    0};
    int length = 3;  
   
    for(int k=0;k< length;k++){
    
        /*R(K-1,K)*/
        if(k>=1){
    
            for(int i=0;i<k;i++){
    
                for(int j=0;j<length;j++){
    
                    R[i][k] += Q[j][i]*init[j][k];
                }
            }
        }
        /*V*/
        for(int j=0;j<length;j++){
    
            if(k<1){
    
             v[j] = init[j][0];
            } else {
    
                    v[j] = init[j][k];
                    for(int g=0;g<k;g++){
    
                        v[j] -= R[g][k]*Q[j][g];
                    }
            }
        }

        
        printf("Jcy v：\n");
        printf("v=\n");
        for(int i=0;i<SIZE;++i)
        {
    
           
                printf("%2.4f ",v[i]);
          
        }
        /*R(K,K)*/
        for(int i=0;i<length;i++){
    
         R[k][k] += v[i]*v[i];
        }
        R[k][k] = sqrt(R[k][k]);
        /*Q**/
        for(int i=0;i<length;i++){
    
            Q[i][k] = v[i]/R[k][k];
        }  


        printf("Jcy QR分解：\n");
        printf("Q=\n");
    for(int i=0;i<SIZE;++i)
    {
    
        for(int j=0;j<SIZE;++j)
        {
    
            printf("%2.4f ",Q[i][j]);
        // 
        }
        printf("\n");
    }

    printf("Jcy QR分解：\n");
    printf("R=\n");
    for(int i=0;i<SIZE;++i)
    {
    
        for(int j=0;j<SIZE;++j)
        {
    
            printf("%2.4f ",R[i][j]);
        // 
        }
        printf("\n");
    }




       //return;
          //Q[k-1][k] = 
    }
    /* printf("Jcy QR分解：\n"); printf("Q=\n"); for(int i=0;i<SIZE;++i) { for(int j=0;j<SIZE;++j) { printf("%2.4f ",Q[i][j]); // } printf("\n"); } printf("Jcy QR分解：\n"); printf("R=\n"); for(int i=0;i<SIZE;++i) { for(int j=0;j<SIZE;++j) { printf("%2.4f ",R[i][j]); // } printf("\n"); } */




}