2015-05-05 57 views
1

所以我有兩個矩陣;一個是2×2,另一種是通過2 1。我想使用的外部庫LAPACK求解線性系統,它好像我需要調用函數dgesv_(),這是使用Lapack解決C中的矩陣

Ax = B 

我必須爲x解決。

所以我真的很困惑他們說的調用函數以及如何將它翻譯成我現在擁有的兩個數組。

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


static double Angstroms[2]; 
static double Energy[2]; 
static double ax[2][2]; 

void file_input(); 
void polynomial(); 


int main() { 

    file_input(); 
    polynomial(); 
    return 0; 
} 

void file_input() { 

    float a, b; 
    int i; 

    FILE * in_file = fopen("H2Mini.txt", "r"); 
    if (outfile == NULL) { 
     printf ("Error file does not exist"); 
     exit (-1); 
    } 
    for (i = 0; i <= 1; i++) { 
      fscanf(in_file, "%f %f\n", &a, &b); 
      Angstroms[i] = a; 
      Energy [i] = b; 
    } 
    fclose(in_file); 

} 

void polynomial() { 

     int i; 
     FILE * outfile = fopen("PolyTest2.txt", "w"); 
     if (outfile == NULL) { 
     printf ("Error file does not exist"); 
     exit (-1); 
     } 
    for (i = 0; i <= 1; i++) { 
     ax[i][0] = 1; 
     fprintf (outfile, "%.8f ", ax[i][0]); 

    } 
    fprintf (outfile, "\n"); 
    for (i = 0; i <= 1; i ++) { 
     ax[i][1] = Angstroms[i]; 
     fprintf (outfile, "%.8f ", ax[i][1]); 
    } 
    } 

所以,我有我的文件是這樣的

[2.00000000 3.00000000 
    6.00000000 5.00000000] 

斧子陣列是要這個樣子

[1.00000000 1.00000000 
    2.00000000 6.00000000] 

和能量陣列是這樣

[3.00000000 5.00000000] 

我的ax數組是Ax = b中的Ax項方程和我的b項是能量陣列。

我看過他們的函數文檔,它在實現它時有點混亂。 在dgesv (n, nrhs, a, lda, ipiv, b, ldb, info)

這個代碼的任何真正明確的例子會超級有用!

+0

「所以我做了兩個矩陣;一種是2 ** 2 **,另一個是** 1乘1 **「。 - 你想要這些? – user2357112

+0

對不起,我的意思是求解矩陣得到x的值爲 –

+0

哪個是矩陣A,哪個x和哪個B? – ztik

回答

0

我不確定'斧頭陣列看起來像這樣'是什麼意思? Ax應該是一個向量不應該嗎?

但是一般調用dgesv使用C包裝LAPACKE應該是這樣的:

LAPACKE_dgesv(
      LAPACK_ROW_MAJOR, // The storage format you use, usually row major in c. 
      n,    // Dimensions of A (nxn) 
      1,    // Number of b's, here always one because you want to solve for a single right hand side. 
      A,    // The input matrix A, it will be destroyed in the call so best make a copy. 
      n,    // LDA basically the length of each column of A in memory. Normally this is also n if A is nxn. 
      p,    // An additional array where lapacke saves pivot ordering. Probably not of any interest to you. Just give it an array of size n. 
      bToX,    // On input input the rhs b, on output this will be the solution vector x 
      1);