CBLAS/LAPACK을 사용하여 C의 대칭 행렬 반전

Question

행렬 및 벡터 곱셈이 필요한 알고리즘을 C로 작성했습니다. I는 매트릭스를 가지고 Q (W X W) 행렬 자체 갖는 벡터 J (1 x 폭)의 전치를 승산과 가산함으로써 생성되는은 I가 사용하여 스케일링 스칼라 .

Q = [(J^T) * J + aI].

I는 벡터 M 얻을 벡터 G와 Q의 역수를 곱해야한다.

M = (Q^(- 1)) * G.

내 알고리즘을 개발 cblas 및 clapack를 사용하고 있습니다. 행렬 Q 는 난수 (float 형)을 사용하여 채워 루틴을 sgetrf_ 및 sgetri_를 사용하여 반전 될 때, 계산 된 역수 올바른이다.

그러나 행렬 Q가 (J^T) x J를 곱할 때 대칭 인 인 경우 계산 된 역함은 잘못되었습니다!.^C에서 LAPACK 루틴을 호출하는 동안

는 I 배열의 행 주요 (C)에서와 열 주요 (FORTRAN에서) 포맷을 알고 있지만, 대칭 행렬이이 같은 문제가 아니다 T = A.

아래의 행렬 반전을위한 C 코드를 첨부했습니다.

나는 이것을 해결하는 더 좋은 방법이있을 것이라고 확신한다. 아무도 이것으로 나를 도울 수 있습니까?

cblas를 사용하는 솔루션이 좋은 것 ...

감사합니다.

void InverseMatrix_R(float *Matrix, int W) 
{ 
    int  LDA = W; 
    int  IPIV[W]; 
    int  ERR_INFO; 
    int  LWORK = W * W; 
    float Workspace[LWORK]; 

    // - Compute the LU factorization of a M by N matrix A 
    sgetrf_(&W, &W, Matrix, &LDA, IPIV, &ERR_INFO); 

    // - Generate inverse of the matrix given its LU decompsotion 
    sgetri_(&W, Matrix, &LDA, IPIV, Workspace, &LWORK, &ERR_INFO); 

    // - Display the Inverted matrix 
    PrintMatrix(Matrix, W, W); 

} 

void PrintMatrix(float* Matrix, int row, int colm) 
{ 
    int i,k; 

    for (i =0; i < row; i++) 
    { 
     for (k = 0; k < colm; k++) 
     { 
      printf("%g, ",Matrix[i*colm + k]); 
     } 

     printf("\n"); 
    } 
} 

Answer 1

나는 BLAS 나 LAPACK을 모른다. 그래서 나는이 행동을 일으키는 원인이 될지 모른다.

그러나 주어진 형식의 행렬에 대해서 역행렬을 계산하는 것은 매우 쉽습니다. 이에 대한 중요한 사실은 구성 요소가 실제 경우 <u|v>이 (내적을 의미

(J^T*J)^2 = (J^T*J)*(J^T*J) = J^T*(J*J^T)*J = <J|J> * (J^T*J) 

입니다 - 정규 선형에게 복잡한 구성 요소 양식을,하지만 당신은 아마 전치하지만 켤레 전치을하지 생각 하는데요 , 그리고 당신은 내면의 제품으로 되돌아 갈 것입니다).

일반화,

(J^T*J)^n = (<J|J>)^(n-1) * (J^T*J), for n >= 1. 

은 우리가 S 및 q에 의해 스칼라 <J|J>에 의해 대칭 정방 행렬 (J^T*J)을 표시하자. 그 후, 충분히 큰 절대 값 ( |a| > q)의 a != 0 일반용 : 사용

(a*I + S) * (1/a*I - 1/(a*(a+q))*S) = I + 1/a*S - 1/(a+q)*S - 1/(a*(a+q))*S^2 
            = I + 1/a*S - 1/(a+q)*S - q/(a*(a+q))*S 
            = I + ((a+q) - a - q)/(a*(a+q))*S 
            = I 

을 계산하여 검증 할 수있는 수식, a = 0 및 a = -q 제외한 모든 a 위해 (analyticity에서) 보관 유지

(a*I + S)^(-1) = 1/a * (I + a^(-1)*S)^(-1) 
       = 1/a * (I + ∑ (-1)^k * a^(-k) * S^k) 
          k>0 
       = 1/a * (I + (∑ (-1)^k * a^(-k) * q^(k-1)) * S) 
          k>0 
       = 1/a * (I - 1/(a+q)*S) 
       = 1/a*I - 1/(a*(a+q))*S 

S^2 = q*S.

또한 계산은 LU 분해를 찾는 것보다 훨씬 간단하고 효율적입니다. 3 × 3 행렬 반전을위한

Answer 2

예, 대한 sgetri.f를 방문 더

//__CLPK_integer is typedef of int 
//__CLPK_real is typedef of float 


__CLPK_integer ipiv[3]; 
{ 
    //Compute LU lower upper factorization of matrix 
    __CLPK_integer m=3; 
    __CLPK_integer n=3; 
    __CLPK_real *a=(float *)this->m1; 
    __CLPK_integer lda=3; 
    __CLPK_integer info; 
    sgetrf_(&m, &n, a, &lda, ipiv, &info); 
} 

{ 
    //compute inverse of a matrix 
    __CLPK_integer n=3; 
    __CLPK_real *a=(float *)this->m1; 
    __CLPK_integer lda=3; 
    __CLPK_real work[3]; 
    __CLPK_integer lwork=3; 
    __CLPK_integer info; 
    sgetri_(&n, a, &lda, ipiv, work, &lwork, &info); 
} 

Answer 3

당신은 LAPACK에 대한 C++ 래퍼를 쉽게 사용할 수있는 Armadillo을 시도 할 수 있습니다. 대칭적인 양의 정부 호 매트릭스

pinv(), 역행렬

solve(), 선형 방정식의 시스템을 해결 (즉 될 수 위해

• inv() 일반적 반전, 선택적 고속화와 함께 : 여러 역 관련 기능을 제공한다 과소 또는 과소 결정), 실제 역관계를하지 않고