C语言-追赶法解方程组（LU分解）

#include <stdio.h> #include <stdlib.h> #include <math.h> #define SUCCESS 1 #define FAIL 0 typedef struct Compoent { struct Compoent *pHead; struct Compoent *pTail; float value; int row; int column; } CompoentNode; CompoentNode *NodeQuery(CompoentNode *head,int row,int column); int DeleteLink(CompoentNode *head); CompoentNode *LU(CompoentNode *head); CompoentNode *RunAfter(CompoentNode* LUHead,CompoentNode *AHead,CompoentNode *b); int main() { printf("*************************************************\n"); printf("* 追赶法求解AX=b V1.1 *\n"); printf("* 于晨 *\n"); printf("* 2015-10-14 *\n"); printf("* 程序可自动识别矩阵维数 *\n"); printf("* 但不会对矩阵A是否对角占优进行检测 *\n"); printf("* 每行以回车结束 *\n"); printf("*************************************************\n"); printf("* 请输入向量矩阵A *\n"); printf("*************************************************\n"); while(1) { printf("A=\n"); //初始化头节点 CompoentNode *head=NULL; CompoentNode *current=NULL; head=(CompoentNode*)malloc(sizeof(CompoentNode)); head->value=0; head->pHead=NULL; head->pTail=NULL; head->column=0; head->row=0; current=head; if(head==NULL) { printf("There is an error: head=(CompoentNode*)malloc(sizeof(CompoentNode)); "); } int row_index=0; int column_index=0; int column_max=0; //接收用户输入矩阵 while(row_index!=column_max||column_max==0) { do { CompoentNode *temp=(CompoentNode*)malloc(sizeof(CompoentNode)); scanf("%f",&(temp->value)); temp->pTail=NULL; temp->pHead=current; temp->row=row_index; temp->column=column_index; current->pTail=temp; current=temp; // printf("%f\n",current->value); head->value++; // printf("%d\n",head->column); column_index++; if(column_max>0) { if(column_index>column_max) break; } }while(getchar()!='\n'); if(row_index==0) { column_max=column_index; } row_index++; column_index=0; head->row=row_index; head->column=column_max; } printf("*************************************************\n"); printf("* 请输入向量b *\n"); printf("*************************************************\n"); printf("b=\n"); CompoentNode *bHead=(CompoentNode*)malloc(sizeof(CompoentNode)); bHead->column=1; bHead->row=head->row; bHead->pHead=NULL; bHead->pTail=NULL; CompoentNode *c=bHead; for(int r=0;r<bHead->row;r++) { CompoentNode *node=(CompoentNode*)malloc(sizeof(CompoentNode)); scanf("%f",&(node->value)); c->pTail=node; node->column=0; node->row=r; node->pTail=NULL; c=node; } printf("*************************************************\n"); printf("LU=\n"); CompoentNode *LUNode=LU(head); for(int r=0;r<head->row;r++) { for(int c=0;c<head->column;c++) { printf("%f ",NodeQuery(LUNode,r,c)->value); } printf("\n"); } printf("*************************************************\n"); printf("x=\n"); CompoentNode *xhead=RunAfter(LUNode,head,bHead); CompoentNode *x=xhead->pTail; while(x!=NULL) { printf("%f\n",x->value); x=x->pTail; } DeleteLink(head); printf("\n"); getchar(); printf("*********************************************************\n"); printf("*********************是否继续? Y/N***********************\n"); if(getchar()=='N') break; } return 0; } /** \根据行列坐标，检索链表中相应节点。找到返回节点指针，未找到返回空 * * \head 头节点 * \row 行坐标 * \column 列坐标 * \return 节点指针 * */ CompoentNode *NodeQuery(CompoentNode *head,int row,int column) { CompoentNode *node=head->pTail; while(node!=NULL) { if(node->row==row&&node->column==column) return node; node=node->pTail; } return NULL; } int DeleteLink(CompoentNode *head) { if(head==NULL) return FAIL; while(head->pTail != NULL) { CompoentNode * p = head->pTail; head->pTail = head->pTail->pTail; free(p); } free(head); return SUCCESS; } /** \brief 根据输入矩阵链表进行LU分解 * * \param head:输入矩阵链表头指针 * \return 返回LU分解后矩阵 * */ CompoentNode *LU(CompoentNode *head) { CompoentNode *current=NULL; CompoentNode *h=(CompoentNode*)malloc(sizeof(CompoentNode)); h->pTail=NULL; h->pHead=NULL; h->value=head->value; h->row=head->row; h->column=head->column; // CompoentNode *temp=head; // while(temp->pTail!=null) current=h; for(int i=0;i<h->row;i++) { if(i==0) { for(int c=0;c<h->column;c++) { CompoentNode *temp=(CompoentNode*)malloc(sizeof(CompoentNode)); temp->value=NodeQuery(head,0,c)->value; temp->column=c; temp->row=0; current->pTail=temp; current=temp; // printf(" %f",LU[c]); } for(int r=1;r<h->row;r++) { // printf(" %f",LU[r*(column_max+1)]); float a=NodeQuery(head,0,0)->value; float value=NodeQuery(head,r,0)->value; CompoentNode *temp=(CompoentNode*)malloc(sizeof(CompoentNode)); temp->value=value/a; temp->column=0; temp->row=r; current->pTail=temp; current=temp; } } else { //处理行 for(int c=i;c<h->column;c++) { float sum=0; for(int k=0;k<=i-1;k++) { float u=NodeQuery(h,k,c)->value; float l=NodeQuery(h,i,k)->value; sum+=l*u; } // printf(" %f",sum); float a=NodeQuery(head,i,c)->value; CompoentNode *temp=(CompoentNode*)malloc(sizeof(CompoentNode)); temp->value=a-sum; temp->column=c; temp->row=i; current->pTail=temp; current=temp; } //处理列 for(int r=i+1;r< h->row;r++) { float sum=0; for(int k=0;k<=i-1;k++) { float u=NodeQuery(h,k,i)->value; float l=NodeQuery(h,r,k)->value; sum+=l*u; } // printf(" %f",sum); float a=NodeQuery(head,r,i)->value; float u=NodeQuery(h,i,i)->value; CompoentNode *temp=(CompoentNode*)malloc(sizeof(CompoentNode)); temp->value=(a-sum)/u; temp->column=i; temp->row=r; current->pTail=temp; current=temp; } } } return h; } //追赶法解方程组 CompoentNode *RunAfter(CompoentNode* LUHead,CompoentNode *AHead,CompoentNode *b) { CompoentNode *yHead=(CompoentNode*)malloc(sizeof(CompoentNode)); yHead->pTail=NULL; yHead->row=b->row; yHead->column=b->column;