Project1

标题: 【C语言发布】非齐次线性方程组的解……（上） [打印本页]

---

作者: zaiy2863    时间: 2016-9-16 18:26
标题: 【C语言发布】非齐次线性方程组的解……（上）


在@RyanBern 的指导下随手写了一发、、、……
但是最高只能支持到阶数n=51 再继续增加会有误差的说(๑• . •๑)
编译器gcc

基本思路是采用高斯消元法、的说
……测试函数保留着，但是已经注释掉了

以及、最喜欢@⑨姐姐 了……
希望⑨姐姐不要嫌弃叶子太笨……(๑• . •๑)

代码如下：

C 代码复制

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
void test_mat_gen(int n, double *a, double *b);
int main(){
        //输入待计算的矩阵  
        int n;
        scanf("%d",&n);
         //n = 84; 
        double A[n*n];//第j列第i行写作A[j*n + i] 
        int i,j,k;
        double b[n+1]; 
        for(i = 0; i<n*n; i++){
                A[i] = 0.0;
        }
        for(i = 0; i<n+1; i++){
                b[i] = 0.0;
        }
//test_mat_gen(n, A, b);//测试函数 
        //录入
        for(i = 0; i<n*n; i++){
                scanf("%lf", &A[i]);
        }
        for(i = 0; i<n; i++){
                scanf("%lf", &b[i]);
        }
        //高斯消元法：对于每个i,j 
        double q;
        k = 0; 
        for(k = 0; k<n; k++){//从A[k*n +k]开始计数 
                for(j = k+1; j<n; j++){//列增加，以j计数 
                q =(double) A[j*n + k]/ A[k*n + k] ;
                //if(q)printf("%f\n", q);
                        for(i = 0; i<n; i++){//j为列 i为行 则固定行数，对列数增加 
                                A[j*n + i] -= q * A[k*n + i];
                        }
                        b[j] -= q* b[k];
                }
        }
 
        //得出结果
         for(i = n-1; i>=0; i--){//对于每个b[i]，减去已经算得的b[i]，用于存储x[i]值 
                 for(j = i+1; j>i; j--){
                 b[i] -= b[j]*A[i*n + j];//b[i]写在第i行         
                 }
                 b[i] /= A[i * n + i];
         }
 
         for(i = 0; i<n; i++){
                printf("%f\n", b[i]); 
        }
 
return 0;        
} 
/* 
void test_mat_gen(int n, double *a, double *b){
        int i;
        for(i = 0; i<n; i++){
                a[i * n + i] = 6;
                if(i < n - 1){
                        a[(i+1)*n + i] = 1;
                }
                if(i){
                        a[(i-1)*n + i] = 8;
                }
                if(i == 0){
                        b[i] = 14;
                }else if(i == n-1){
                        b[i] = 7;
                }else{
                        b[i] = 15;
                }
        }
 
}
*/

#include <stdio.h>


#include <math.h>


#include <stdlib.h>


void test_mat_gen(int n, double *a, double *b);


int main(){


        //输入待计算的矩阵  


        int n;


        scanf("%d",&n);


         //n = 84; 


        double A[n*n];//第j列第i行写作A[j*n + i] 


        int i,j,k;


        double b[n+1]; 


        for(i = 0; i<n*n; i++){


                A[i] = 0.0;


        }


        for(i = 0; i<n+1; i++){


                b[i] = 0.0;


        }


//test_mat_gen(n, A, b);//测试函数 


        //录入


        for(i = 0; i<n*n; i++){


                scanf("%lf", &A[i]);


        }


        for(i = 0; i<n; i++){


                scanf("%lf", &b[i]);


        }


        //高斯消元法：对于每个i,j 


        double q;


        k = 0; 


        for(k = 0; k<n; k++){//从A[k*n +k]开始计数 


                for(j = k+1; j<n; j++){//列增加，以j计数 


                q =(double) A[j*n + k]/ A[k*n + k] ;


                //if(q)printf("%f\n", q);


                        for(i = 0; i<n; i++){//j为列 i为行 则固定行数，对列数增加 


                                A[j*n + i] -= q * A[k*n + i];


                        }


                        b[j] -= q* b[k];


                }


        }


 


        //得出结果


         for(i = n-1; i>=0; i--){//对于每个b[i]，减去已经算得的b[i]，用于存储x[i]值 


                 for(j = i+1; j>i; j--){


                 b[i] -= b[j]*A[i*n + j];//b[i]写在第i行         


                 }


                 b[i] /= A[i * n + i];


         }


 


         for(i = 0; i<n; i++){


                printf("%f\n", b[i]); 


        }


 


return 0;        


} 


/* 


void test_mat_gen(int n, double *a, double *b){


        int i;


        for(i = 0; i<n; i++){


                a[i * n + i] = 6;


                if(i < n - 1){


                        a[(i+1)*n + i] = 1;


                }


                if(i){


                        a[(i-1)*n + i] = 8;


                }


                if(i == 0){


                        b[i] = 14;


                }else if(i == n-1){


                        b[i] = 7;


                }else{


                        b[i] = 15;


                }


        }


 


}


*/

---

作者: RyanBern    时间: 2016-9-16 19:01
前排。坐等更新。

另外这代码还得注意程序规范的问题，比如声明的时候通常不加变量名字，只说类型。
还有这个写法到底是啥意思

C 代码复制

scanf("%.0f\t", &[i]);

scanf("%.0f\t", &[i]);

---

作者: 寒冷魔王    时间: 2016-9-17 01:17


哇，小叶子也做这个了呢~
前排支持！

另外，看你使用一维数组作为矩阵处理的话，可以考虑写一些对应的函数哦~

1. double A[n*n];//第j列第i行写作A[j*n + i]
   
2. // 写成：
   
3. double getMatrixElement(double matrix[], int i, int j) // 按照习惯交换j和i的位置，根据你的注释，j表示列，i表示行，数学上的习惯是a(i j)
   
4. {
   
5.     return matrix[j*n + i];
   
6. }
   
7. void setMatrixElement(double matrix[], int i, int j, double element)
   
8. {
   
9.     matrix[j*n + i] = element;
   
10. }
    

复制代码

用法：
double A[n * n];
v = getMatrixElement(A, i, j);
setMatrixElement(A, i, j, v);
这样写不容易出错，只要把i和j的位置对应好就行了~


============================以下内容可以无视==============================

小叶子你的写法呢，double A[n*n]; 其实是C99的标准，这个n是变量。
如果要用C++(兼容C89但不兼容C99)编译呢，最好写成

1. double *A = (double*)calloc(n*n, sizeof(double));  // 用calloc的话不用手动清零。

复制代码

当然需要用完后free才行。

建议使用struct来表示Matrix

1. struct Matrix
   
2. {
   
3.     double *data;
   
4. };

复制代码

用struct Matrix代替double[]的好处是不会和其他的double[]弄混，而且能够直接复制struct Matrix哦~

1. void initMatrix(struct Matrix *matrix, int n)
   
2. {
   
3.     matrix->data = (double*)calloc(n*n, sizeof(double));
   
4. }
   
5. void freeMatrix(struct Matrix *matrix)
   
6. {
   
7.     free(matrix->data);
   
8. }
   
9. double getMatrixElement(struct Matrix *matrix, int i, int j)
   
10. {
    
11.     return matrix->data[j*n + i];
    
12. }
    
13. void setMatrixElement(struct Matrix *matrix, int i, int j, double element)
    
14. {
    
15.     matrix->data[j*n + i] = element;
    
16. }
    

复制代码

使用方法是
struct Matrix A;
initMatrix(&A, n);              // 初始化
v = getMatrixElement(&A, i, j); // 获取 v = A(i j)
setMatrixElement(&A, i, j, v);  // 设置 A(i j) = v
freeMatrix(&A);                 // 释放


如果在struct Matrix内加个变量储存矩阵大小，还可以顺便检查其访问元素是否越界。

1. struct Matrix
   
2. {
   
3.     double *data;
   
4.     int size;
   
5. };

复制代码

更好一些的代码：

1. void initMatrix(struct Matrix *matrix, int n)
   
2. {
   
3.     matrix->data = (double*)calloc(n*n, sizeof(double));
   
4.     matrix->size = n;
   
5. }
   
6. void freeMatrix(struct Matrix *matrix)
   
7. {
   
8.     free(matrix->data);
   
9. }
   
10. double getMatrixElement(struct Matrix *matrix, int i, int j)
    
11. {
    
12.     assert(i >= 0 && matrix->size > i && j >= 0 && matrix->size > j);
    
13.     return matrix->data[j*n + i];
    
14. }
    
15. void setMatrixElement(struct Matrix *matrix, int i, int j, double element)
    
16. {
    
17.     assert(i >= 0 && matrix->size > i && j >= 0 && matrix->size > j);
    
18.     matrix->data[j*n + i] = element;
    
19. }
    

复制代码

其中，getMatrixElement和setMatrixElement新增的assert会在Debug模式下检查不符合表达式的错误，自动结束程序。
它检查的是 i 和 j 在 [0,matrix->size) 的区间内。
assert用于检查编程错误，比如循环时不小心多了1次，越界了之类的。不要用assert判断程序需要处理的错误。

这个assert在非Debug模式下会去除，GCC中，使用NDEBUG用于去除assert。
因为assert会去除，所以不要在assert中写能够改变程序状态的代码，如
assert(i++);
这会导致程序的行为不一致。

---

作者: RyanBern    时间: 2016-9-18 00:01

寒冷魔王 发表于 2016-9-17 01:17
哇，小叶子也做这个了呢~
前排支持！

首先，是因为是叶子所以回复这么认真么/w\

如果要用C++(兼容C89但不兼容C99)编译呢，最好写成...

其实我早就跟叶子说了我不建议这样，但是谁让它是一个 C 语言程序并且自动打开 gcc 的 c99 标准呢……似乎我好像永远活在了 ANSI C 标准中无法自拔。

另外“使用C++编译”似乎是个不准确的说法呢……

用struct Matrix代替double[]的好处是不会和其他的double[]弄混，而且能够直接复制struct Matrix哦~

你确定？除非你显式写出复制的函数，否则是绝对不能这样“直接”复制的。

如果单纯考虑不弄混的话，可以直接定义一个类型，这样就不容易混了。

C 代码复制

typedef double* MatElePr;

typedef double* MatElePr;

另外从你的回帖里我能看出你尽量把程序做得面向对象一点，于是想方设法去写 Constructor Destructor 等东西。其实我在群里是这样教叶子的：不许用 struct，不许有面向对象的思维，给我用纯指针弄。

看来我的反面向对象思维越来越严重了。不过，有关的高性能的数学库确实不会考虑这些封装呢。

---

作者: 寒冷魔王    时间: 2016-9-18 00:50

RyanBern 发表于 2016-9-18 00:01
首先，是因为是叶子所以回复这么认真么/w\

/w\感谢您的回复，看到小叶子这么认真所以激动的有些忘乎所以了呢~

关于

你确定？除非你显式写出复制的函数，否则是绝对不能这样“直接”复制的。

我是确定的，刚刚测试了一下（gcc -std=c99）

看来我没记错，是可以直接复制的。（应该是新标准？好像是当年看C++Primer里面写的，具体出处我也忘了）

不能弄混是指传递参数的时候啦，用typedef定义的话传递参数还是会弄混的（当然R君要是不用函数的话确实能够达到效果）

另外呢，可能是我经常用C++的缘故，所以用C的话有些面向对象了呢（/w\）
哎嘿嘿，正统的面向过程C语言确实不喜欢面向对象什么的。。
其实因为我用到了calloc和free，这种构建使用函数来管理是更方便及安全的。（比如只是init而忘了free，导致内存泄漏，也可以通过在函数内部加输出判断是否泄漏）
同时，关于您对性能的顾虑，我觉得您多虑了。
现代编译器（比如小叶子手头上的gcc（目测5.1以上））的话，能够正确处理这些问题，自动inline什么的，编译出来的效果可以媲美C语言的宏。
而且inline函数可以进行参数检查，可以预防很多编程错误。要比直接写效果更好。
哦对，inline貌似是C99标准。。。

我觉得初次编写的时候，这个安全性要首先保证。这样编写和修改都很简单。
编写成型以后再做优化（比如内嵌汇编神马的），这样可以避免过早优化。

其实我写的时候就想，要是直接用C++的RAII多方便。
如果使用C++的话，我觉得基于RAII的内存管理和适当的封装对于编写清晰安全并且高效的代码很有帮助。像小叶子写的A[j*n + i]，改成A.get(i, j)什么的，我觉得可读性就上升许多，也易于维护。不需要手动调用free什么的，也减少了内存泄漏的可能性。

所以说适当的封装是有好处的，R君可能需要常年编写这种高性能的程序，对这方面比较敏感。不过我觉得这种损失目前来说还是可以忽略的。

总之！
像R君这样的数学专业人士、人品保障的前辈，能给出专业意见、提出人生建议，呵护小叶子健康快乐成长，我对您是很放心的。

---

作者: 晴兰    时间: 2016-9-18 01:46
提示: 作者被禁止或删除 内容自动屏蔽

---

作者: RyanBern    时间: 2016-9-18 08:43

寒冷魔王 发表于 2016-9-18 00:50
/w\感谢您的回复，看到小叶子这么认真所以激动的有些忘乎所以了呢~

关于

C 代码复制

struct A1{
  double data[5];
};

struct A1{


  double data[5];


};

C 代码复制

struct A2{
  double *data;
};

struct A2{


  double *data;


};

下面两个有本质的区别啊，不要把它们当成一样的了喂！

C 代码复制

printf("%d %d\n", sizeof(struct A1), sizeof(struct A2));

printf("%d %d\n", sizeof(struct A1), sizeof(struct A2));

这样会看到很明显的不同。

概况来说，如果是动态分配的内存，是绝对不可以这样玩的。

C 代码复制

struct A{
  double *data;
};
 
int main(int argc, char **argv){
  struct A a1, a2;
  /* suppose we call 'constructors' for a1 and a2 */
  init(&a1, n); init(&a2, n);
  a2 = a1; /* the problem occurs ! */
  a1.data[0] = 1; /* then what is a2.data[0]? */
  /* suppose we call 'destructors' for a1 and a2 */
  finalize(&a1); finalize(&a2); /* core dumped !*/
  return 0;
}

struct A{


  double *data;


};


 


int main(int argc, char **argv){


  struct A a1, a2;


  /* suppose we call 'constructors' for a1 and a2 */


  init(&a1, n); init(&a2, n);


  a2 = a1; /* the problem occurs ! */


  a1.data[0] = 1; /* then what is a2.data[0]? */


  /* suppose we call 'destructors' for a1 and a2 */


  finalize(&a1); finalize(&a2); /* core dumped !*/


  return 0;


}

写得有点仓促，但是我觉得意思应该表示清楚了。

---

作者: condir    时间: 2016-9-18 10:40
我看了半天也没有想通，为什么在RPGMaker发布这么专业的C语言技术文档
能问一句，这个是做什么用的？用在哪里？难道RPGMaker里还能用C语言吗？

---

作者: 喵呜喵5    时间: 2016-9-18 13:32

condir 发表于 2016-9-18 10:40
我看了半天也没有想通，为什么在RPGMaker发布这么专业的C语言技术文档
能问一句，这个是做什么用的 ...

因为真正能用在 RPG Maker 里的东西都发布到水区去了

---

作者: 浮云半仙    时间: 2016-11-17 14:43


膝行而前，莫能仰视

---

作者: RyanBern    时间: 2016-11-17 19:31
糊了一个 PLU，里面只有 PLU 分解，没有解方程组的部分。

叶子看完了可以把解方程组的部分补上。不知道有没有 BUG，没仔细看过。

源码：
https://github.com/RyanBernX/math/tree/master/lu
说明书：
https://ryanbernx.github.io/math/

---

欢迎光临 Project1 (https://rpg.blue/)

Powered by Discuz! X3.1