2015-05-05 53 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);