Q&A
「3.i行j列の成分が0ならば、n[i]=0,非0ならばn[i]=1を格納。これをj列の成分すべてについて行う。」これを文面通りに受け取ると、最終列の結果がn[]に格納されていることになるんですが、おそらくそういうことじゃないと思います。
###前提・実現したいこと
Cで行列のランク計算をするプログラムを作っています。
プログラムの流れは、
1.行列の行数(lnum)、列数(cnum)を入力
2.行列の各成分を入力、格納。
3．行列の走査（matsearch）を行う。
走査範囲は初めは1列目の全行であるとする。
１列目から順に
特定の列の走査範囲内の要素をを上から順に確認
→すべて0 ⇒ 走査範囲はそのままで次の列に
→上端が非0かつ上端以外はすべて0 ⇒ 走査範囲を一つ縮める（元の走査範囲から上端の行を抜いたもの）。そして走査範囲は次の列へ
→上端以外に非0が存在
⇒ その列が基本変形対象範囲最重要列J₀
かつ
その走査範囲の上端の行は基本変形対象範囲開始行I₀。
for行ループ → break
for列ループ → break
matsearch → 終了
最終的に…matsearchの返却値は
①J=cnumですべてその列の成分が０
→ -2
②（I₀，J₀）＝０かつbreak
→ I₀
③（I₀，J₀）≠０かつbreak
→ -1
4.行基本変形
①：
rank計算へ
②：
（I, J₀）成分のうち、I₀≦I≦lnumの範囲で
非0となっている成分について考える。
このような成分のうち、最小行番号をもつ成分の行番号
をI₁とする。（I₀＜I₁≦lnum）
2-2-1.1：（I₁，J）成分でI₁行成分をすべて割る。
2-2-1.2：（I₀，J）成分に
「（I₀，J）成分 + （I₁，J）成分」
の演算結果を代入する。
2-2-1.3：
（I，J）成分に（I＞I₀）
「（I，J）成分 -｛（I，J）成分*（I₀，J）｝」
の演算結果を代入する。
そのまま③へ
③：
操作2-1
2-1.1：（I₀ , J₀）成分でI₀行成分をすべて割る。
2-1.2：（I，J）成分に（I＞I₀）
「（I，J）成分 -｛（I，J）成分*（I₀，J）｝」
の演算結果を代入する。
という手順を踏んで行基本変形を行う予定だったのですが、この手順通りに操作がいかないです。
どうかよろしくおねがいいいたします。
###発生している問題・エラーメッセージ
計算過程でループが終了せずどこが間違っているのかがわからないです。
###該当のソースコード
#include <stdio.h> #define NA 7 #define NB 3 #define PLUS 0 #define MINUS 1 #define ERROR_R 1.000000e-15 int zerojudgef(double a) { return ((a < ERROR_R) & (a > (-1)*ERROR_R) ? 0 : 1); } int matsearchf(double mat[][NA][NA], int lnum, int cnum, int keylinenum, int keycolumnnum, int rank) /* keylinenum : the number of line that is the start point of elementary row operation. Exactly, in the column of key-column number, and below the line of the start line of elemntary row operation ,there may exists the element that is not zero which has the smallest number of line. The line number of such elements is defined as keylinenum. That is, if the start line of elementary row operation is not the exactly defined keyline number, the start line of elementray row operation is linebreak. keycolumnnum : the number of column that is the key-column of elementary row op ertion. */ { int i, j; int linebreak = 0; for(j=0; j<cnum; j++){ int k = 0; int notzero = 0; int keynum_for_count = 0; for(i=linebreak; i<lnum; i++){ if(zerojudgef(mat[0][i][j])){ notzero++; }/*for-for-if*/ }/*for-for*/ /*Finding out the key-line number In the column of key-column number, and below the line of key-line number, there may exists the element that is not zero which has the smallest number of line. The line number of such elements is defined as keylinenum, and in this function, for the purpose of counting the key-line number, a variable 'keynum_for_count' is defined. */ do{ keynum_for_count = (zerojudgef(mat[0][k][j]) != 0 ? k : -1);//ここは最悪k+linebreakでもok k++; }while( (keynum_for_count < 0) & (k != (NA-1) ) );//matのイランとこは全部0で初期化が条件 if(zerojudgef(notzero)){/*There exist at least one element that is not equal to zero.*/ if(zerojudgef(mat[0][linebreak][j])){/*The top elements in the scope of line-searching is not equal to zero.*/ if(notzero > 1){/*There exists not-zero elemnet at the top of the scope of line searching, and exists not-zero elements below the elements at the top. break roop*/ keylinenum = keynum_for_count; keycolumnnum = j; return -1; break; }else{ /*There exists not-zero element only at the top of the searching range, and elements below it is all equal to zero. continue roop*/ linebreak++; }/*for-if-if-if*/ }else{/*The top elements in the scope of line-searching is equal to zero, but there exist at least one elements that is not equal to zero. break roop*/ keylinenum = keynum_for_count; keycolumnnum = j; return linebreak; break; }/*for-if-if*/ }else{/*There exist only elements that is equal to zero. continue roop*/ }/*for-if*/ }/*for*/ rank = linebreak; return -2; /*Finally, the result is classified into 3 types. 1:All the elemnets in the last column of the scope is equal to zero. output = -2 2:The (keylinenum, keycolumnnum) elemnet is equal to zero. output = linebreak 3:The (keylinenum, keycolumnnum) elemnet is not equal to zero. output = -1 */ } void divf(int divdo, int divdonel, int cnum, double mat[][NA][NA]) { /* divdo :割る数の要素の列番号 divdonel :割られる数の要素が存在する行番号 cnum :割り算実行される行の列数 */ double divisor = mat[0][divdonel][divdo]; int i; for(i=0; i<cnum; i++){ mat[0][divdonel][i] /= divisor; } /**/ int line, column; puts("divf"); printf("\n"); for(line=0; line<NA; line++){ printf("|"); for(column=0; column<NA; column++){ printf("%5.2f", mat[0][line][column]); } printf(" |\n"); } printf("\n\n"); /**/ } void subtraf(int p, int q, int r, int cnum, double mat[][NA][NA], double coeff) { /* This function operates the line p ± the line q. Judge by whether r = 0 or 1 -r＝０→addition＋ -r＝１→subtractionー cnum = The total number of line coeff = the coefficient of the number that subtacts any numbers */ int pm = (r == 1 ? (-1) : 1); int i; for(i=0; i<cnum; i++){ mat[0][p][i] += (coeff * pm * mat[0][q][i]); } /**/ int line, column; puts("subtraf"); printf("\n"); for(line=0; line<NA; line++){ printf("|"); for(column=0; column<NA; column++){ printf("%5.2f", mat[0][line][column]); } printf(" |\n"); } printf("\n\n"); /**/ } void erowoperatef(double mat[][NA][NA], int lnum, int cnum, int keylinenum, int keycolumnnum) { int rank = 0; int i, j, k; int breakjudge = 0; for(i=0; i<NA; i++){ k = ( ( matsearchf(mat, lnum, cnum, keylinenum, keycolumnnum, rank) < 0 ) ? matsearchf(mat, lnum, cnum, keylinenum, keycolumnnum, rank) : 1 ); switch(k){ case 1:{ divf(keycolumnnum, keylinenum, cnum, mat); subtraf( matsearchf(mat, lnum, cnum, keylinenum, keycolumnnum, rank), keylinenum, PLUS, cnum, mat, 1); }/*for-switch-case0*/ case -1:{ divf(keycolumnnum, keylinenum, cnum, mat); if(i = keylinenum){ }else{ for(j=0; j<lnum; j++){ subtraf( j, keylinenum, MINUS, cnum, mat , mat[0][j][keycolumnnum]); }/*for-switch-case-1-if-for*/ }/*for-switch-case-1-if*/ break; }/*for-switch-case1*/ case -2:{ breakjudge = 1; break; }/*for-switch-case-1*/ }/*switch*/ if(breakjudge == 1){ break; }/*for-if*/ }/*for*/ } int rankcaloperatef(double mat[][NA][NA], int lnum, int cnum) { int keylinenum = 0; int keycolumnnum = 0; erowoperatef(mat, lnum, cnum, keylinenum, keycolumnnum); int rank = 0; int empty = matsearchf(mat, lnum, cnum, keylinenum, keycolumnnum, rank); return rank; } int main(void) { double mat[NB][NA][NA] = {}; int lnum, cnum; /*mat[kind of matrix][line number][column number] */ puts("\nWe will calculate the rank of a matrix.\nPlease enter the scale of matrix."); printf("The number of lines (under %d) = ", NA); scanf("%d", &lnum); printf("The number of columns (under %d) = ", NA); scanf("%d", &cnum); puts("Please enter the elements of the matrix."); puts("'mat[i][j]' represents the (i,j) elements of the matrix."); int i, j; for(i=0; i<lnum; i++){ for(j=0; j<cnum; j++){ printf("mat[%d][%d] = ", i+1, j+1); scanf("%lf", &mat[0][i][j]); mat[1][i][j] = mat[0][i][j]; } } /*mat[1]にはもともと行列が保存される。*/ int Rank = rankcaloperatef(mat, lnum, cnum); int r; for(r=0; r<2; r++){ printf("\n"); for(i=0; i<lnum; i++){ printf("|"); for(j=0; j<cnum; j++){ printf("%5.f", mat[r][i][j]); } printf(" |\n"); } printf("\n\n"); } printf("\nThe rank of this matrix is %d.\n", Rank); return 0; }
###試したこと
上のソースコードのなかで、２つの/**/で挟まれたprintf …の部分は演算過程を見るために付け加えたものです。いくつかのデバッグは自己解決できたのですが、これに関してはわからないです。
又、初めのアルゴリズムでは欠陥があったので、少し変えて作り直しました。
例えば「100. j=0..Nについて手順101~103を繰り返す」みたいに正確に書いてもらえますか?
回答2件
0
自己解決
#問題点
*matsearchf内
erowoperatefやrankcaloperatefでmatsearchfを使用する時に渡す仮引数にkeylinenum,keycolumnnum,rankがあった。matsearchfでこれらの仮引数に数値を代入して、それをmatsearchf外のブロック有効範囲でも利用する予定だった。しかし、ポインタ型変数に指定なったために、これらの仮引数に代入した数値はmatsearchfのプロック有効範囲しかもたず、matsearchf外では当然のことながらこれらの変数が元の関数内で保持していた値に戻ってしまう。その為、正しくない値で以降の計算を行ってしまっていた。
→これらの3変数をポインタ型変数として扱うことにより解決。
*erowoperatef内
C
1 case -1:{ 2 divf(keycolumnnum, keylinenum, cnum, mat); 3 if(i = keylinenum){ 4 }else{ 5 for(j=0; j<lnum; j++){ 6 subtraf( j, keylinenum, MINUS, cnum, mat , mat[0][j][keycolumnnum]); 7 }/*for-switch-case-1-if-for*/ 8 }/*for-switch-case-1-if*/ 9 break; 10 }/*for-switch-case1*/ 11
の部分で「行ー行」を行う際に「I₀行ーI₀行」を行わないようにとif文で「引かれる側の行」の行番号を判断できるようにしていた。しかし、
-for文で行番号として使用している変数がｊであるのに対し、変数iをif文での判定で用いていて無意味なことを行っていた。
→書き換え。
-if文の判定式で等価演算子==を用いていた筈であったが、=が１つだけで代入操作を行っていた。
→書き換え。
-for文とif文の前後の位置関係を間違えていた。
本来は
C
1 2 for(j=0; j<lnum; j++){ 3 if(j == keylinenum){ 4 … 5 } 6 } 7
とすべきだった。
→書き替え。
#完成版のソースコード
C
1#include <stdio.h> 2#define NA 7 3#define NB 3 4#define PLUS 0 5#define MINUS 1 6#define ERROR_R 1.000000e-15 7 8 9int zerojudgef(double a) 10{ 11 return ((a < ERROR_R) & (a > (-1)*ERROR_R) ? 0 : 1); 12} 13 14void divf(int divdo, int divdonel, int cnum, double mat[][NA][NA]) 15{ 16 /* 17 divdo :割る数の要素の列番号 18 divdonel :割られる数の要素が存在する行番号 19 cnum :割り算実行される行の列数 20 */ 21 double divisor = mat[0][divdonel][divdo]; 22 23 int i; 24 for(i=0; i<cnum; i++){ 25 mat[0][divdonel][i] /= divisor; 26 } 27 28 29} 30 31void subtraf(int p, int q, int r, int cnum, double mat[][NA][NA], double coeff) 32{ 33 /* 34 This function operates the line p ± the line q. 35 Judge by whether r = 0 or 1 36 -r＝０→addition＋ 37 -r＝１→subtractionー 38 cnum = The total number of line 39 coeff = the coefficient of the number that subtacts any numbers 40 */ 41 int pm = (r == 1 ? (-1) : 1); 42 int i; 43 for(i=0; i<cnum; i++){ 44 mat[0][p][i] += (coeff * pm * mat[0][q][i]); 45 } 46 47 48 49} 50 51 52int matsearchf(double mat[][NA][NA], int lnum, int cnum, int * keylinenum, int * keycolumnnum, int * rank) 53/* 54keylinenum : the number of line that is the start point of elementary row operation. Exactly, in the column of key-column number, and below the line of the start line of elemntary row operation ,there may exists the element that is not zero which has the smallest number of line. The line number of such elements is defined as keylinenum. 55That is, if the start line of elementary row operation is not the exactly defined keyline number, the start line of elementray row operation is linebreak. 56 57keycolumnnum : the number of column that is the key-column of elementary row op 58ertion. 59*/ 60{ 61 int i, j; 62 int linebreak = 0; 63 for(j=0; j<cnum; j++){ 64 int k = 0; 65 int notzero = 0; 66 int keynum_for_count = 0; 67 for(i=linebreak; i<lnum; i++){ 68 69 if(zerojudgef(mat[0][i][j])){ 70 notzero++; 71 }/*for-for-if*/ 72 }/*for-for*/ 73/*Finding out the key-line number 74In the column of key-column number, and below the line of key-line number, there may exists the element that is not zero which has the smallest number of line. The line number of such elements is defined as keylinenum, and in this function, for the purpose of counting the key-line number, a variable 'keynum_for_count' is defined. 75*/ 76 do{ 77 keynum_for_count = (zerojudgef(mat[0][k+linebreak][j]) != 0 ? (k+linebreak) : -1);//ここは最悪k+linebreakでもok 78/* 79If 'zerojudgef(mat[0][k][j])' is equal to zero, 'keynum_for_count' (= -1) must not be substituted into variable 'keylinenum'. This is because serching range will continue to be same as before. 80 81If 'zerojudgef(mat[0][k][j])' is not equal to zero, 'keynum_for_count' must be substituted into variable 'keylinenum'. 82*/ 83 k++; 84 }while( (keynum_for_count < 0) & ((k+linebreak) != (NA-1) ) );//matのイランとこは全部0で初期化が条件 85 if(zerojudgef(notzero)){ 86/*There exist at least one element that is not equal to zero.*/ 87 88 if(zerojudgef(mat[0][linebreak][j])){ 89 /*The top elements in the scope of line-searching is not equal to zero.*/ 90 if(notzero > 1){ 91 /*There exists not-zero elemnet at the top of the scope of line searching, and exists not-zero elements below the elements at the top. break roop*/ 92 *keylinenum = keynum_for_count; 93 *keycolumnnum = j; 94 *rank = linebreak; 95 return -1; 96 break; 97 }else{ 98/*There exists not-zero element only at the top of the searching range, and elements below it is all equal to zero. continue roop*/ 99 100 linebreak++; 101 *keylinenum = keynum_for_count; 102 *keycolumnnum = j; 103 *rank = linebreak; 104 105 106 }/*for-if-if-if*/ 107 }else{ 108 /*The top elements in the scope of line-searching is equal to zero, but there exist at least one elements that is not equal to zero. break roop*/ 109 *keylinenum = keynum_for_count; 110 *keycolumnnum = j; 111 *rank = linebreak; 112 return linebreak; 113 break; 114 }/*for-if-if*/ 115 116 }else{ 117/*There exist only elements that is equal to zero. continue roop*/ 118 }/*for-if*/ 119 }/*for*/ 120 *rank = linebreak; 121 return -2; 122 123 /* 124Finally, the result is classified into 3 types. 1251:All the elemnets in the last column of the scope is equal to zero. 126output = -2 127for example, 128|_ 1_ 0_ 0__| 129|_ 0_ 1_ 0__| 130|_ 0_ 0_ 0__| 131 1322:The (keylinenum, keycolumnnum) elemnet is equal to zero. 133output = linebreak 1343:The (keylinenum, keycolumnnum) elemnet is not equal to zero. 135output = -1 136*/ 137} 138 139void erowoperatef(double mat[][NA][NA], int lnum, int cnum) 140{ 141 int keylinenum = 0; 142 int keycolumnnum = 0; 143 int rank = 0; 144 int i, j, k; 145 int breakjudge = 0; 146 for(i=0; i<NA; i++){ 147 k = ( ( matsearchf(mat, lnum, cnum, &keylinenum, &keycolumnnum, &rank) < 0 ) ? matsearchf(mat, lnum, cnum, &keylinenum, &keycolumnnum, &rank) : 1 ); 148 switch(k){ 149 case 1:{ 150 divf(keycolumnnum, keylinenum, cnum, mat); 151 subtraf( matsearchf(mat, lnum, cnum, &keylinenum, &keycolumnnum, &rank), keylinenum, PLUS, cnum, mat, 1); 152 153 }/*for-switch-case0*/ 154 case -1:{ 155 156 divf(keycolumnnum, keylinenum, cnum, mat); 157 for(j=0; j<lnum; j++){ 158 if(j == keylinenum){ 159 }else{ 160 subtraf( j, keylinenum, MINUS, cnum, mat , mat[0][j][keycolumnnum]); 161 }/*for-switch-case-1-for-if*/ 162 }/*for-switch-case-1-for*/ 163 164 break; 165 }/*for-switch-case1*/ 166 case -2:{ 167 /*Go to calcuration of rank.*/ 168 169 /*Finding out the element that has the smallest column number in lowest line that has the not-zero element */ 170 int topdivcol = 0; 171 do{ 172 if(zerojudgef(mat[0][(rank-1)][topdivcol])){ 173 break; 174 }else{ 175 topdivcol++; 176 } 177 }while(topdivcol < cnum); 178 179 /*Operating the elementary row operation to the lowest line that has not-zero elements, and topdivcolumn.*/ 180 if(topdivcol == cnum){ 181 }else{ 182 183 divf(topdivcol, (rank-1), cnum, mat); 184 for(j=0; j<lnum; j++){ 185 if(j == (rank-1)){ 186 }else{ 187 subtraf( j, (rank-1), MINUS, cnum, mat , mat[0][j][(rank-1)]); 188 }/*for-switch-case-2-if-for-if*/ 189 }/*for-switch-case-2-if-for*/ 190 191 } 192 193 breakjudge = 1; 194 break; 195 }/*for-switch-case-1*/ 196 }/*switch*/ 197 if(breakjudge == 1){ 198 199 200 break; 201 }/*for-if*/ 202 }/*for*/ 203 204} 205 206 207 208 209 210 211int rankcaloperatef(double mat[][NA][NA], int lnum, int cnum) 212{ 213 int keylinenum = 0; 214 int keycolumnnum = 0; 215 erowoperatef(mat, lnum, cnum); 216 int Rank = 0; 217 int empty = matsearchf(mat, lnum, cnum, &keylinenum, &keycolumnnum, &Rank); 218 return Rank; 219} 220 221 222void matprintf(double mat[][NA][NA], int matnumber, int lnum, int cnum) 223{ 224 225int i, j; 226 227if(matnumber == 1){printf("Original Matrix\n"); 228}else{ 229printf("Operated Matrix\n"); 230} 231 232 233for(i=0; i<lnum; i++){ 234 printf("|_"); 235 for(j=0; j<cnum; j++){ 236 printf("%5.2f_", mat[matnumber][i][j]); 237 } 238 printf("_|\n\n"); 239 } 240 241} 242 243 244int main(void) 245{ 246 247 static double mat[NB][NA][NA] = {}; 248 int lnum, cnum; 249 /*mat[kind of matrix][line number][column number] */ 250 puts("\nWe will calculate the rank of a matrix.\nPlease enter the scale of matrix."); 251 printf("The number of lines (under %d) = ", NA); 252 scanf("%d", &lnum); 253 printf("The number of columns (under %d) = ", NA); 254 scanf("%d", &cnum); 255 puts("Please enter the elements of the matrix."); 256 puts("'mat[i][j]' represents the (i,j) elements of the matrix."); 257 int i, j; 258 for(i=0; i<lnum; i++){ 259 for(j=0; j<cnum; j++){ 260 printf("mat[%d][%d] = ", i+1, j+1); 261 scanf("%lf", &mat[0][i][j]); 262 mat[1][i][j] = mat[0][i][j]; 263 } 264 } 265 /*Original matrix is input into mat[1].*/ 266 267 int Rank = rankcaloperatef(mat, lnum, cnum); 268 matprintf(mat, 1, lnum, cnum); 269 matprintf(mat, 0, lnum, cnum); 270 printf("\nThe rank of this matrix is %d.\n", Rank); 271 return 0; 272} 273 274
趣味でCを勉強し始めたばかりの初心者でございます故、初歩的なミスがかなり多かったです。
お付き合いいただきありがとうございました。
投稿2017/02/16 00:47
総合スコア16
0
計算誤差により "i行j列の成分が0ならば" の判定が狂うことはありませんか?
一般に float/doubleなどの整数型でないものを ==, != で判定してはならないのですが。
投稿2017/02/10 07:00
総合スコア16614
回答ありがとうございます。
ご指摘の通り、浮動小数点型を==で判定していました。
そのため、浮動小数点型を判定できるように修正を行いました。
しかし、根本的には解決できなかったため、もう一度書き直しましたが上手くいっておりません。
あなたの回答
tips
プレビュー
