c语言矩阵运算库（c++矩阵运算库函数）

1.一个优雅好用的c语言库必须满足哪些条件
2.实现一个矩阵运算库的几点思考
- (1)采用预定义的数据类型，避免直接使用编译器定义的数据类型
- (2)基于对象编程，定义矩阵对象
- (3)除了特别编写的内存处理函数(使用栈链表保存、释放动态分配的内存地址)，不允许任何函数直接分配和释放内存
- (4)防御性编程，对输入参数做有效性检查，并返回错误号
- (5)注意编程细节的打磨
3.完整c程序
参考资料

编程既是技术输出也是艺术创作。鉴赏高手写的程序，往往让人眼前一亮，他们思路、逻辑清晰，所呈现的代码简洁、优雅、高效，令人为之叹服。烂代码则像“屎山"一般让人作呕，软件难以维护最大的原因除了需求模糊的客观因素，最重要的主观因素还是代码写得烂。有生之年，愿能保持对编程的热情，不断提升编程能力，真正体会其乐趣，共勉！

1.一个优雅好用的c语言库必须满足哪些条件

这里给出的软件开发应遵循的一般原则，摘自Les Piegl和Wayne Tiller所著的一本非常经典的书The NURBS Book(Second Edition)。
(1)工具性(toolability):应该使用可用的工具来建立新的应用程序；
(2)可移植性(portability):应用程序应该容易被移植到不同的软件和硬件平台；
(3)可重用性(reusability):程序的编写应该便于重复使用代码段；
(4)可检验性(testability):程序应该简单一致，以便于测试和调试；
(5)可靠性(reliability):对程序运行过程中可能出现的各种错误应该进行合理、一致的处理，以使系统稳定、可靠；
(6)可扩展性(enhanceability):代码必须易于理解，以便可以容易地增加新的功能；
(7)可维护性(fixability):易于找出程序错误的位置；
(8)一致性(consistency):在整个库中，编程习惯应保持一致；
(9)可读性(communicability):程序应该易于阅读和理解;
(10)编程风格(style of programming):代码看起来像书上的数学公式那样以便于读者理解，同时遵循用户友好的编程风格;
(11)易用性(usability):应该使一些非专业的用户也能够方便地使用所开发的库来开发各种更高层次的应用程序;
(12)数值高效性(numerical efficiency):所有函数必须仔细推敲，保证其数值高效性;
(13)基于对象编程(object based programming):避免在函数间传递大量数据，并且使代码易于理解。

2.实现一个矩阵运算库的几点思考

(1)采用预定义的数据类型，避免直接使用编译器定义的数据类型

1 2 3 4 5 6 7 typedef unsigned int ERROR_ID; typedef int INDEX; typedef short FLAG; typedef int INTEGER; typedef double REAL; typedef char* STRING; typedef void VOID;

使用预定义的数据类型，有利于程序移植，且可以提高可读性。例如，如果一个系统只支持单精度浮点数，那么只需修改数据类型REAL为float，达到一劳永逸的效果。定义INDEX与INTEGER数据类型是为了在编程时区分索引变量与普通整形变量，同样提高可读性。

(2)基于对象编程，定义矩阵对象

1 2 3 4 5 6 typedef struct matrix { INTEGER rows; INTEGER columns; REAL* p; }MATRIX;

这里，用一级指针而非二级指针指向矩阵的数据内存地址，有诸多原因，详见博文:为什么我推荐使用一级指针创建二维数组？。

(3)除了特别编写的内存处理函数(使用栈链表保存、释放动态分配的内存地址)，不允许任何函数直接分配和释放内存

不恰当的分配、使用与释放内存很可能导致内存泄漏、系统崩溃等致命的错误。如果一个函数需动态申请多个内存，那么可能会写出这样啰嗦的程序:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 double* x = NULL, * y = NULL, * z = NULL; x = (double*)malloc(n1 * sizeof(double)); if (x == NULL) return -1; y = (double*)malloc(n2 * sizeof(double)); if (y == NULL) { free(x); x = NULL; return -1; } z = (double*)malloc(n3 * sizeof(double)); if (z == NULL) { free(x); x = NULL; free(y); y = NULL; return -1; }

为了优雅地实现动态内存分配与释放，Les Piegl大神分3步来处理内存申请与释放：
a)在进入一个新的程序时，一个内存堆栈被初始化为空；
b)当需要申请内存时，调用特定的函数来分配所需的内存，并将指向内存的指针存入堆栈中的正确位置；
c)在离开程序时，遍历内存堆栈，释放其中的指针所指向的内存。
程序结构大致如下：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 STACKS S; MATRIX* m = NULL; INTEGER rows = 3, columns = 4; ERROR_ID errorID = _ERROR_NO_ERROR; init_stack(&S); m = creat_matrix(rows, columns, &errorID, &S); if (m == NULL) goto EXIT; //do something // ... EXIT: free_stack(&S); return errorID;

(4)防御性编程，对输入参数做有效性检查，并返回错误号

例如输入的矩阵行数、列数应该是正整数，指针参数必须非空等等。

(5)注意编程细节的打磨

a)操作符(逗号，等号等)两边必须空一格；
b)逻辑功能相同的程序间不加空行，逻辑功能独立的程序间加空行；
c)条件判断关键字(for if while等)后必须加一空格，起到强调作用，也更清晰；
d)函数内部定义局部变量后，必须空一行后再编写函数主体。

3.完整c程序

本矩阵运算库只包含了矩阵的基本运算，包括创建任意二维/三维矩阵、创建零矩阵及单位矩阵、矩阵加法、矩阵减法、矩阵乘法、矩阵求逆、矩阵转置、矩阵的迹、矩阵LUP分解、解矩阵方程AX=B。

common.h

matrix.h

common.c

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 /******************************************************************************* * File Name : common.c * Library/Module Name : MatrixComputation * Author : Marc Pony(marc_pony@163.com) * Create Date : 2021/7/16 * Abstract Description : 矩阵运算库公用源文件 *******************************************************************************/ /******************************************************************************* * (1)Debug Switch Section *******************************************************************************/ /******************************************************************************* * (2)Include File Section *******************************************************************************/ #include "common.h" /******************************************************************************* * (3)Macro Define Section *******************************************************************************/ /******************************************************************************* * (4)Struct(Data Types) Define Section *******************************************************************************/ /******************************************************************************* * (5)Prototype Declare Section *******************************************************************************/ /******************************************************************************* * (6)Global Variable Declare Section *******************************************************************************/ /******************************************************************************* * (7)File Static Variable Define Section *******************************************************************************/ /******************************************************************************* * (8)Function Define Section *******************************************************************************/ /********************************************************************************************** Function: init_stack Description: 初始化栈 Input: 无 Output: 无 Input_Output: 栈指针 Return: 无 Author: Marc Pony(marc_pony@163.com) ***********************************************************************************************/ VOID init_stack(_IN_OUT STACKS* S) { if (S == NULL) { return; } memset(S, 0, sizeof(STACKS)); } /********************************************************************************************** Function: free_stack Description: 释放栈 Input: 栈指针 Output: 无 Input_Output: 无 Return: 无 Author: Marc Pony(marc_pony@163.com) ***********************************************************************************************/ VOID free_stack(_IN STACKS* S) { MATRIX_NODE* matrixNode = NULL; MATRIX_ELEMENT_NODE* matrixElementNode = NULL; if (S == NULL) { return; } while (S->matrixNode != NULL) { matrixNode = S->matrixNode; S->matrixNode = matrixNode->next; free(matrixNode->ptr); matrixNode->ptr = NULL; free(matrixNode); matrixNode = NULL; } while (S->matrixElementNode != NULL) { matrixElementNode = S->matrixElementNode; S->matrixElementNode = matrixElementNode->next; free(matrixElementNode->ptr); matrixElementNode->ptr = NULL; free(matrixElementNode); matrixElementNode = NULL; } // ... // 释放其他指针 } matrix.c /******************************************************************************* * File Name : matrix.c * Library/Module Name : MatrixComputation * Author : Marc Pony(marc_pony@163.com) * Create Date : 2021/2/24 * Abstract Description : 矩阵运算库源文件 *******************************************************************************/ /******************************************************************************* * (1)Debug Switch Section *******************************************************************************/ /******************************************************************************* * (2)Include File Section *******************************************************************************/ #include "matrix.h" /******************************************************************************* * (3)Macro Define Section *******************************************************************************/ /******************************************************************************* * (4)Struct(Data Types) Define Section *******************************************************************************/ /******************************************************************************* * (5)Prototype Declare Section *******************************************************************************/ /******************************************************************************* * (6)Global Variable Declare Section *******************************************************************************/ /******************************************************************************* * (7)File Static Variable Define Section *******************************************************************************/ /******************************************************************************* * (8)Function Define Section *******************************************************************************/ VOID print_matrix(MATRIX* a, STRING string) { INDEX i, j; printf("matrix %s:", string); printf("\n"); for (i = 0; i < a->rows; i++) { for (j = 0; j < a->columns; j++) { printf("%f ", a->p[i * a->columns + j]); } printf("\n"); } printf("\n"); } /********************************************************************************************** Function: creat_matrix Description: 创建矩阵 Input: 矩阵行数rows，列数columns Output: 错误号指针errorID，栈指针S Input_Output: 无 Return: 矩阵指针 Author: Marc Pony(marc_pony@163.com) ***********************************************************************************************/ MATRIX* creat_matrix(_IN INTEGER rows, _IN INTEGER columns, _OUT ERROR_ID* errorID, _OUT STACKS* S) { MATRIX* matrix = NULL; MATRIX_NODE* matrixNode = NULL; MATRIX_ELEMENT_NODE* matrixElementNode = NULL; *errorID = _ERROR_NO_ERROR; if (rows <= 0 || columns <= 0 || errorID == NULL || S == NULL) { *errorID = _ERROR_INPUT_PARAMETERS_ERROR; return NULL; } matrix = (MATRIX*)malloc(sizeof(MATRIX)); matrixNode = (MATRIX_NODE*)malloc(sizeof(MATRIX_NODE)); matrixElementNode = (MATRIX_ELEMENT_NODE*)malloc(sizeof(MATRIX_ELEMENT_NODE)); if (matrix == NULL || matrixNode == NULL || matrixElementNode == NULL) { free(matrix); matrix = NULL; free(matrixNode); matrixNode = NULL; free(matrixElementNode); matrixElementNode = NULL; *errorID = _ERROR_FAILED_TO_ALLOCATE_HEAP_MEMORY; return NULL; } matrix->rows = rows; matrix->columns = columns; matrix->p = (REAL*)malloc(rows * columns * sizeof(REAL)); //确保matrix非空才能执行指针操作 if (matrix->p == NULL) { free(matrix->p); matrix->p = NULL; free(matrix); matrix = NULL; free(matrixNode); matrixNode = NULL; free(matrixElementNode); matrixElementNode = NULL; *errorID = _ERROR_FAILED_TO_ALLOCATE_HEAP_MEMORY; return NULL; } matrixNode->ptr = matrix; matrixNode->next = S->matrixNode; S->matrixNode = matrixNode; matrixElementNode->ptr = matrix->p; matrixElementNode->next = S->matrixElementNode; S->matrixElementNode = matrixElementNode; return matrix; } /********************************************************************************************** Function: creat_multiple_matrices Description: 创建多个矩阵 Input: 矩阵行数rows，列数columns，个数count Output: 错误号指针errorID，栈指针S Input_Output: 无 Return: 矩阵指针 Author: Marc Pony(marc_pony@163.com) ***********************************************************************************************/ MATRIX* creat_multiple_matrices(_IN INTEGER rows, _IN INTEGER columns, _IN INTEGER count, _OUT ERROR_ID* errorID, _OUT STACKS* S) { INDEX i; MATRIX* matrix = NULL, *p = NULL; MATRIX_NODE* matrixNode = NULL; *errorID = _ERROR_NO_ERROR; if (rows <= 0 || columns <= 0 || count <= 0 || errorID == NULL || S == NULL) { *errorID = _ERROR_INPUT_PARAMETERS_ERROR; return NULL; } matrix = (MATRIX*)malloc(count * sizeof(MATRIX)); matrixNode = (MATRIX_NODE*)malloc(sizeof(MATRIX_NODE)); if (matrix == NULL || matrixNode == NULL) { free(matrix); matrix = NULL; free(matrixNode); matrixNode = NULL; *errorID = _ERROR_FAILED_TO_ALLOCATE_HEAP_MEMORY; return NULL; } for (i = 0; i < count; i++) { p = creat_matrix(rows, columns, errorID, S); if (p == NULL) { free(matrix); matrix = NULL; free(matrixNode); matrixNode = NULL; *errorID = _ERROR_FAILED_TO_ALLOCATE_HEAP_MEMORY; return NULL; } matrix[i] = *p; } matrixNode->ptr = matrix; matrixNode->next = S->matrixNode; S->matrixNode = matrixNode; return matrix; } /********************************************************************************************** Function: creat_zero_matrix Description: 创建零矩阵 Input: 矩阵行数rows，列数columns Output: 错误号指针errorID，栈指针S Input_Output: 无 Return: 矩阵指针 Author: Marc Pony(marc_pony@163.com) ***********************************************************************************************/ MATRIX* creat_zero_matrix(_IN INTEGER rows, _IN INTEGER columns, _OUT ERROR_ID* errorID, _OUT STACKS* S) { MATRIX* matrix = NULL; *errorID = _ERROR_NO_ERROR; if (rows <= 0 || columns <= 0 || errorID == NULL || S == NULL) { *errorID = _ERROR_INPUT_PARAMETERS_ERROR; return NULL; } matrix = creat_matrix(rows, columns, errorID, S); if (*errorID == _ERROR_NO_ERROR) { memset(matrix->p, 0, rows * columns * sizeof(REAL)); } return matrix; } /********************************************************************************************** Function: creat_eye_matrix Description: 创建单位矩阵 Input: 矩阵行数rows，列数columns Output: 错误号指针errorID，栈指针S Input_Output: 无 Return: 矩阵指针 Author: Marc Pony(marc_pony@163.com) ***********************************************************************************************/ MATRIX* creat_eye_matrix(_IN INTEGER n, _OUT ERROR_ID* errorID, _OUT STACKS* S) { INDEX i; MATRIX* matrix = NULL; *errorID = _ERROR_NO_ERROR; if (n <= 0 || errorID == NULL || S == NULL) { *errorID = _ERROR_INPUT_PARAMETERS_ERROR; return NULL; } matrix = creat_matrix(n, n, errorID, S); if (*errorID == _ERROR_NO_ERROR) { memset(matrix->p, 0, n * n * sizeof(REAL)); for (i = 0; i < n; i++) { matrix->p[i * n + i] = 1.0; } } return matrix; } /********************************************************************************************** Function: matrix_add Description: 矩阵A + 矩阵B = 矩阵C Input: 矩阵A,矩阵B Output: 矩阵C Input_Output: 无 Return: 错误号 Author: Marc Pony(marc_pony@163.com) ***********************************************************************************************/ ERROR_ID matrix_add(_IN MATRIX* A, _IN MATRIX* B, _OUT MATRIX* C) { INDEX i, j; INTEGER rows, columns; ERROR_ID errorID = _ERROR_NO_ERROR; if (A == NULL || B == NULL || C == NULL) { errorID = _ERROR_INPUT_PARAMETERS_ERROR; return errorID; } if (A->rows != B->rows || A->rows != C->rows || B->rows != C->rows || A->columns != B->columns || A->columns != C->columns || B->columns != C->columns) { errorID = _ERROR_MATRIX_ROWS_OR_COLUMNS_NOT_EQUAL; return errorID; } rows = A->rows; columns = A->columns; for (i = 0; i < rows; i++) { for (j = 0; j < columns; j++) { C->p[i * columns + j] = A->p[i * columns + j] + B->p[i * columns + j]; } } return errorID; } /********************************************************************************************** Function: matrix_subtraction Description: 矩阵A - 矩阵B = 矩阵C Input: 矩阵A,矩阵B Output: 矩阵C Input_Output: 无 Return: 错误号 Author: Marc Pony(marc_pony@163.com) ***********************************************************************************************/ ERROR_ID matrix_subtraction(_IN MATRIX* A, _IN MATRIX* B, _OUT MATRIX* C) { INDEX i, j; INTEGER rows, columns; ERROR_ID errorID = _ERROR_NO_ERROR; if (A == NULL || B == NULL || C == NULL) { errorID = _ERROR_INPUT_PARAMETERS_ERROR; return errorID; } if (A->rows != B->rows || A->rows != C->rows || B->rows != C->rows || A->columns != B->columns || A->columns != C->columns || B->columns != C->columns) { errorID = _ERROR_MATRIX_ROWS_OR_COLUMNS_NOT_EQUAL; return errorID; } rows = A->rows; columns = A->columns; for (i = 0; i < rows; i++) { for (j = 0; j < columns; j++) { C->p[i * columns + j] = A->p[i * columns + j] - B->p[i * columns + j]; } } return errorID; } /********************************************************************************************** Function: matrix_multiplication Description: 矩阵乘法C = A * B Input: 矩阵A,矩阵B Output: 矩阵C Input_Output: 无 Return: 错误号 Author: Marc Pony(marc_pony@163.com) ***********************************************************************************************/ ERROR_ID matrix_multiplication(_IN MATRIX* A, _IN MATRIX* B, _OUT MATRIX* C) { INDEX i, j, k; REAL sum; ERROR_ID errorID = _ERROR_NO_ERROR; if (A == NULL || B == NULL || C == NULL) { errorID = _ERROR_INPUT_PARAMETERS_ERROR; return errorID; } if (A->columns != B->rows || A->rows != C->rows || B->columns != C->columns) { errorID = _ERROR_MATRIX_MULTIPLICATION; return errorID; } for (i = 0; i < A->rows; i++) { for (j = 0; j < B->columns; j++) { sum = 0.0; for (k = 0; k < A->columns; k++) { sum += A->p[i * A->columns + k] * B->p[k * B->columns + j]; } C->p[i * B->columns + j] = sum; } } return errorID; } /********************************************************************************************** Function: matrix_inverse Description: 矩阵求逆 Input: 矩阵A Output: 矩阵A的逆矩阵 Input_Output: 无 Return: 错误号 Author: Marc Pony(marc_pony@163.com) ***********************************************************************************************/ ERROR_ID matrix_inverse(_IN MATRIX* A, _OUT MATRIX* invA) { INDEX i; INTEGER n; MATRIX* ATemp = NULL; ERROR_ID errorID = _ERROR_NO_ERROR; STACKS S; if (A == NULL || invA == NULL) { errorID = _ERROR_INPUT_PARAMETERS_ERROR; return errorID; } if (A->rows != A->columns) { errorID = _ERROR_MATRIX_MUST_BE_SQUARE; return errorID; } init_stack(&S); n = A->rows; ATemp = creat_matrix(n, n, &errorID, &S); if (errorID != _ERROR_NO_ERROR) goto EXIT; memcpy(ATemp->p, A->p, n * n * sizeof(REAL)); memset(invA->p, 0, n * n * sizeof(REAL)); for (i = 0; i < n; i++) { invA->p[i * n + i] = 1.0; } errorID = solve_matrix_equation_by_lup_decomposition(ATemp, invA); EXIT: free_stack(&S); return errorID; } /********************************************************************************************** Function: matrix_transpose Description: 矩阵转置 Input: 矩阵A Output: 矩阵A的转置 Input_Output: 无 Return: 错误号 Author: Marc Pony(marc_pony@163.com) ***********************************************************************************************/ ERROR_ID matrix_transpose(_IN MATRIX* A, _OUT MATRIX* transposeA) { INDEX i, j; ERROR_ID errorID = _ERROR_NO_ERROR; if (A == NULL || transposeA == NULL) { errorID = _ERROR_INPUT_PARAMETERS_ERROR; return errorID; } if (A->rows != transposeA->columns || A->columns != transposeA->rows) { errorID = _ERROR_MATRIX_TRANSPOSE_FAILED; return errorID; } for (i = 0; i < A->rows; i++) { for (j = 0; j < A->columns; j++) { transposeA->p[j * A->rows + i] = A->p[i * A->columns + j]; } } return errorID; } /********************************************************************************************** Function: matrix_trace Description: 矩阵的迹 Input: 矩阵A Output: 矩阵A的迹 Input_Output: 无 Return: 错误号 Author: Marc Pony(marc_pony@163.com) ***********************************************************************************************/ ERROR_ID matrix_trace(_IN MATRIX* A, _OUT REAL *trace) { INDEX i; ERROR_ID errorID = _ERROR_NO_ERROR; if (A == NULL || trace == NULL) { errorID = _ERROR_INPUT_PARAMETERS_ERROR; return errorID; } if (A->rows != A->columns) { errorID = _ERROR_MATRIX_MUST_BE_SQUARE; return errorID; } *trace = 0.0; for (i = 0; i < A->rows; i++) { *trace += A->p[i * A->columns + i]; } return errorID; } /********************************************************************************************** Function: lup_decomposition Description: n行n列矩阵LUP分解PA = L * U Input: n行n列矩阵A Output: n行n列下三角矩阵L，n行n列上三角矩阵U，n行n列置换矩阵P Input_Output: 无 Return: 错误号 Author: Marc Pony(marc_pony@163.com) 参考：https://zhuanlan.zhihu.com/p/84210687 ***********************************************************************************************/ ERROR_ID lup_decomposition(_IN MATRIX* A, _OUT MATRIX* L, _OUT MATRIX* U, _OUT MATRIX* P) { INDEX i, j, k, index, s, t; INTEGER n; REAL maxValue, temp; ERROR_ID errorID = _ERROR_NO_ERROR; if (A == NULL || L == NULL || U == NULL || P == NULL) { errorID = _ERROR_INPUT_PARAMETERS_ERROR; return errorID; } if (A->rows != A->columns) { errorID = _ERROR_MATRIX_MUST_BE_SQUARE; return errorID; } n = A->rows; memcpy(U->p, A->p, n * n * sizeof(REAL)); memset(L->p, 0, n * n * sizeof(REAL)); memset(P->p, 0, n * n * sizeof(REAL)); for (i = 0; i < n; i++) { L->p[i * n + i] = 1.0; P->p[i * n + i] = 1.0; } for (j = 0; j < n - 1; j++) { //Select i(>= j) that maximizes | U(i, j) | index = -1; maxValue = 0.0; for (i = j; i < n; i++) { temp = fabs(U->p[i * n + j]); if (temp > maxValue) { maxValue = temp; index = i; } } if (index == -1) { continue; } //Interchange rows of U : U(j, j : n) < ->U(i, j : n) for (k = j; k < n; k++) { s = j * n + k; t = index * n + k; temp = U->p[s]; U->p[s] = U->p[t]; U->p[t] = temp; } //Interchange rows of L : L(j, 1 : j - 1) < ->L(i, 1 : j - 1) for (k = 0; k < j; k++) { s = j * n + k; t = index * n + k; temp = L->p[s]; L->p[s] = L->p[t]; L->p[t] = temp; } //Interchange rows of P : P(j, 1 : n) < ->P(i, 1 : n) for (k = 0; k < n; k++) { s = j * n + k; t = index * n + k; temp = P->p[s]; P->p[s] = P->p[t]; P->p[t] = temp; } for (i = j + 1; i < n; i++) { s = i * n + j; L->p[s] = U->p[s] / U->p[j * n + j]; for (k = j; k < n; k++) { U->p[i * n + k] -= L->p[s] * U->p[j * n + k]; } } } return errorID; } /********************************************************************************************** Function: solve_matrix_equation_by_lup_decomposition Description: LUP分解解矩阵方程AX=B,其中A为n行n列矩阵，B为n行m列矩阵，X为n行m列待求矩阵(写到矩阵B) Input: n行n列矩阵A Output: 无 Input_Output: n行m列矩阵B(即n行m列待求矩阵X) Return: 错误号 Author: Marc Pony(marc_pony@163.com) ***********************************************************************************************/ ERROR_ID solve_matrix_equation_by_lup_decomposition(_IN MATRIX* A, _IN_OUT MATRIX* B) { INDEX i, j, k, index, s, t; INTEGER n, m; REAL sum, maxValue, temp; MATRIX* L = NULL, * U = NULL, * y = NULL; ERROR_ID errorID = _ERROR_NO_ERROR; STACKS S; if (A == NULL || B == NULL) { errorID = _ERROR_INPUT_PARAMETERS_ERROR; return errorID; } if (A->rows != A->columns) { errorID = _ERROR_MATRIX_MUST_BE_SQUARE; return errorID; } init_stack(&S); n = A->rows; m = B->columns; L = creat_matrix(n, n, &errorID, &S); if (errorID != _ERROR_NO_ERROR) goto EXIT; U = creat_matrix(n, n, &errorID, &S); if (errorID != _ERROR_NO_ERROR) goto EXIT; y = creat_matrix(n, m, &errorID, &S); if (errorID != _ERROR_NO_ERROR) goto EXIT; memcpy(U->p, A->p, n * n * sizeof(REAL)); memset(L->p, 0, n * n * sizeof(REAL)); for (i = 0; i < n; i++) { L->p[i * n + i] = 1.0; } for (j = 0; j < n - 1; j++) { //Select i(>= j) that maximizes | U(i, j) | index = -1; maxValue = 0.0; for (i = j; i < n; i++) { temp = fabs(U->p[i * n + j]); if (temp > maxValue) { maxValue = temp; index = i; } } if (index == -1) { continue; } //Interchange rows of U : U(j, j : n) < ->U(i, j : n) for (k = j; k < n; k++) { s = j * n + k; t = index * n + k; temp = U->p[s]; U->p[s] = U->p[t]; U->p[t] = temp; } //Interchange rows of L : L(j, 1 : j - 1) < ->L(i, 1 : j - 1) for (k = 0; k < j; k++) { s = j * n + k; t = index * n + k; temp = L->p[s]; L->p[s] = L->p[t]; L->p[t] = temp; } //Interchange rows of P : P(j, 1 : n) < ->P(i, 1 : n), C = P * B，等价于对B交换行 for (k = 0; k < m; k++) { s = j * m + k; t = index * m + k; temp = B->p[s]; B->p[s] = B->p[t]; B->p[t] = temp; } for (i = j + 1; i < n; i++) { s = i * n + j; L->p[s] = U->p[s] / U->p[j * n + j]; for (k = j; k < n; k++) { U->p[i * n + k] -= L->p[s] * U->p[j * n + k]; } } } for (i = 0; i < n; i++) { if (fabs(U->p[i * n + i]) < 1.0e-20) { errorID = _ERROR_MATRIX_EQUATION_HAS_NO_SOLUTIONS; goto EXIT; } } //L * y = C for (j = 0; j < m; j++) { for (i = 0; i < n; i++) { sum = 0.0; for (k = 0; k < i; k++) { sum = sum + L->p[i * n + k] * y->p[k * m + j]; } y->p[i * m + j] = B->p[i * m + j] - sum; } } //U * x = y for (j = 0; j < m; j++) { for (i = n - 1; i >= 0; i--) { sum = 0.0; for (k = i + 1; k < n; k++) { sum += U->p[i * n + k] * B->p[k * m + j]; } B->p[i * m + j] = (y->p[i * m + j] - sum) / U->p[i * n + i]; } } EXIT: free_stack(&S); return errorID; }

test_matrix.c

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 #include "matrix.h" void main() { REAL a[3 * 3] = { 1,2,3,6,5,5,8,7,2 }; REAL b[3 * 3] = {1,2,3,6,5,4,3,2,1}; MATRIX *A = NULL, * B = NULL, * C = NULL, * D = NULL, * E = NULL, * Z = NULL, * invA = NULL, *m = NULL; ERROR_ID errorID = _ERROR_NO_ERROR; REAL trace; STACKS S; init_stack(&S); Z = creat_zero_matrix(3, 3, &errorID, &S); print_matrix(Z, "Z"); E = creat_eye_matrix(3, &errorID, &S); print_matrix(E, "E"); A = creat_matrix(3, 3, &errorID, &S); A->p = a; print_matrix(A, "A"); B = creat_matrix(3, 3, &errorID, &S); B->p = b; print_matrix(B, "B"); C = creat_matrix(3, 3, &errorID, &S); D = creat_matrix(3, 3, &errorID, &S); invA = creat_matrix(3, 3, &errorID, &S); errorID = matrix_add(A, B, C); errorID = matrix_subtraction(A, B, C); errorID = matrix_multiplication(A, B, C); errorID = matrix_transpose(A, D); print_matrix(D, "D"); errorID = matrix_trace(A, &trace); errorID = matrix_inverse(A, invA); print_matrix(invA, "invA"); m = creat_multiple_matrices(3, 3, 2, &errorID, &S); m[0].p = a; m[1].p = b; free_stack(&S); }

参考资料

The NURBS Book(Second Edition). Les Piegl,Wayne Tiller

到此这篇关于纯c语言优雅地实现矩阵运算库的方法的文章就介绍到这了,更多相关c语言矩阵运算内容请搜索服务器之家以前的文章或继续浏览下面的相关文章希望大家以后多多支持服务器之家！

原文链接：https://blog.csdn.net/maple_2014/article/details/119720296

如果您对该产品感兴趣，请填写办理（客服微信：xiaoxiongyidong）

关于作者：访客

c语言矩阵运算库（c++矩阵运算库函数）

目录

1.一个优雅好用的c语言库必须满足哪些条件

2.实现一个矩阵运算库的几点思考

(1)采用预定义的数据类型，避免直接使用编译器定义的数据类型

(2)基于对象编程，定义矩阵对象

(3)除了特别编写的内存处理函数(使用栈链表保存、释放动态分配的内存地址)，不允许任何函数直接分配和释放内存

(4)防御性编程，对输入参数做有效性检查，并返回错误号

(5)注意编程细节的打磨

3.完整c程序

参考资料

发表评论取消回复

业务办理

运营商合作

目录

1.一个优雅好用的c语言库必须满足哪些条件

2.实现一个矩阵运算库的几点思考

(1)采用预定义的数据类型，避免直接使用编译器定义的数据类型

(2)基于对象编程，定义矩阵对象

(3)除了特别编写的内存处理函数(使用栈链表保存、释放动态分配的内存地址)，不允许任何函数直接分配和释放内存

(4)防御性编程，对输入参数做有效性检查，并返回错误号

(5)注意编程细节的打磨

3.完整c程序

参考资料

为您推荐：

发表评论 取消回复

业务办理

运营商合作

发表评论取消回复