src/jp/nyatla/nyartoolkit/core/NyARMat.java

   1 /* \r
   2  * PROJECT: NyARToolkit\r
   3  * --------------------------------------------------------------------------------\r
   4  * This work is based on the original ARToolKit developed by\r
   5  *   Hirokazu Kato\r
   6  *   Mark Billinghurst\r
   7  *   HITLab, University of Washington, Seattle\r
   8  * http://www.hitl.washington.edu/artoolkit/\r
   9  *\r
  10  * The NyARToolkit is Java version ARToolkit class library.\r
  11  * Copyright (C)2008 R.Iizuka\r
  12  *\r
  13  * This program is free software; you can redistribute it and/or\r
  14  * modify it under the terms of the GNU General Public License\r
  15  * as published by the Free Software Foundation; either version 2\r
  16  * of the License, or (at your option) any later version.\r
  17  * \r
  18  * This program is distributed in the hope that it will be useful,\r
  19  * but WITHOUT ANY WARRANTY; without even the implied warranty of\r
  20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\r
  21  * GNU General Public License for more details.\r
  22  * \r
  23  * You should have received a copy of the GNU General Public License\r
  24  * along with this framework; if not, write to the Free Software\r
  25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA\r
  26  * \r
  27  * For further information please contact.\r
  28  *      http://nyatla.jp/nyatoolkit/\r
  29  *      <airmail(at)ebony.plala.or.jp>\r
  30  * \r
  31  */\r
  32 package jp.nyatla.nyartoolkit.core;\r
  33 \r
  34 import jp.nyatla.nyartoolkit.NyARException;\r
  35 \r
  36  \r
  37 \r
  38 /**\r
  39  *      ARMat構造体に対応するクラス\r
  40  *      typedef struct {\r
  41  *              double *m;\r
  42  *              int row;\r
  43  *              int clm;\r
  44  *      }ARMat;\r
  45  *\r
  46  */\r
  47 public class NyARMat{\r
  48     /**\r
  49      * 配列サイズと行列サイズは必ずしも一致しないことに注意\r
  50      * 返された配列のサイズを行列の大きさとして使わないこと！\r
  51      * \r
  52      */\r
  53     protected double[][] m;\r
  54     private int clm,row;\r
  55     /**\r
  56      * デフォルトコンストラクタは機能しません。\r
  57      * @throws NyARException\r
  58      */\r
  59     protected NyARMat() throws NyARException\r
  60     {\r
  61         throw new NyARException();\r
  62     }\r
  63     public NyARMat(int i_row,int i_clm)\r
  64     {\r
  65         m=new double[i_row][i_clm];\r
  66         clm=i_clm;\r
  67         row=i_row;\r
  68     }\r
  69     /**\r
  70      * i_row x i_clmサイズの行列を格納できるように行列サイズを変更します。\r
  71      * 実行後、行列の各値は不定になります。\r
  72      * @param i_row\r
  73      * @param i_clm\r
  74      */\r
  75     public void realloc(int i_row,int i_clm)\r
  76     {\r
  77         if(i_row<=this.m.length && i_clm<=this.m[0].length)\r
  78         {\r
  79             //十分な配列があれば何もしない。\r
  80         }else{\r
  81             //不十分なら取り直す。\r
  82             m=new double[i_row][i_clm];\r
  83         }\r
  84         this.clm=i_clm;\r
  85         this.row=i_row;\r
  86     }\r
  87     public int getClm()\r
  88     {\r
  89         return clm;\r
  90     }\r
  91     public int getRow()\r
  92     {\r
  93         return row;\r
  94     }\r
  95     /**\r
  96      * 行列をゼロクリアする。\r
  97      */\r
  98     public void zeroClear()\r
  99     {\r
 100         int i,i2;\r
 101         //For順変更OK\r
 102         for(i=row-1;i>=0;i--){\r
 103             for(i2=clm-1;i2>=0;i2--){\r
 104                 m[i][i2]=0.0;\r
 105             }\r
 106         }\r
 107     }\r
 108     /**\r
 109      * i_copy_fromの内容を自分自身にコピーします。\r
 110      * 高さ・幅は同一で無いと失敗します。\r
 111      * @param i_copy_from\r
 112      */\r
 113     public void copyFrom(NyARMat i_copy_from)throws NyARException\r
 114     {\r
 115         //サイズ確認\r
 116         if(this.row!=i_copy_from.row ||this.clm!=i_copy_from.clm)\r
 117         {\r
 118             throw new NyARException();\r
 119         }\r
 120         //値コピー\r
 121         for(int r=this.row-1;r>=0;r--){\r
 122             for(int c=this.clm-1;c>=0;c--){\r
 123                 this.m[r][c]=i_copy_from.m[r][c];\r
 124             }\r
 125         }\r
 126     }\r
 127 \r
 128     public double[][] getArray()\r
 129     {\r
 130         return m;\r
 131     }\r
 132 //    public void getRowVec(int i_row,NyARVec o_vec)\r
 133 //    {\r
 134 //      o_vec.set(this.m[i_row],this.clm);\r
 135 //    }\r
 136     /**\r
 137      * aとbの積を自分自身に格納する。arMatrixMul()の代替品\r
 138      * @param a\r
 139      * @param b\r
 140      * @throws NyARException\r
 141      */\r
 142     public void matrixMul(NyARMat a, NyARMat b) throws NyARException\r
 143     {\r
 144         if(a.clm != b.row || this.row != a.row || this.clm != b.clm){\r
 145             throw new NyARException();\r
 146         }\r
 147         double w;\r
 148         int r,c,i;\r
 149         double[][] am=a.m,bm=b.m,dm=this.m;\r
 150         //For順変更禁止\r
 151         for(r = 0; r < this.row; r++){\r
 152             for(c = 0; c < this.clm; c++){\r
 153                 w=0.0;//dest.setARELEM0(r, c,0.0);\r
 154                 for(i = 0; i < a.clm; i++){\r
 155                     w+=am[r][i]*bm[i][c];//ARELEM0(dest, r, c) += ARELEM0(a, r, i) * ARELEM0(b, i, c);\r
 156                 }\r
 157                 dm[r][c]=w;\r
 158             }\r
 159         }\r
 160     }\r
 161     private int[] wk_nos_matrixSelfInv=new int[50];\r
 162 //    private final static double matrixSelfInv_epsl=1.0e-10;\r
 163     /**\r
 164      * i_targetを逆行列に変換する。arMatrixSelfInv()と、arMatrixSelfInv_minv()関数を合成してあります。\r
 165      * OPTIMIZE STEP[485->422]\r
 166      * @param i_target\r
 167      * 逆行列にする行列\r
 168      * @return\r
 169      * 逆行列があればTRUE/無ければFALSE\r
 170      * \r
 171      * @throws NyARException\r
 172      */\r
 173     public boolean matrixSelfInv() throws NyARException\r
 174     {\r
 175         double[][] ap=this.m;\r
 176         int dimen=this.row;\r
 177         int dimen_1=dimen-1;\r
 178         double[] ap_n,ap_ip,ap_i;//wap;\r
 179         int j,ip,nwork;\r
 180         int[] nos=wk_nos_matrixSelfInv;//この関数で初期化される。\r
 181         //double epsl;\r
 182         double p,pbuf,work;\r
 183 \r
 184         /* check size */\r
 185         switch(dimen){\r
 186         case 0:\r
 187             throw new NyARException();\r
 188         case 1:\r
 189             ap[0][0]=1.0/ap[0][0];//*ap = 1.0 / (*ap);\r
 190             return true;/* 1 dimension */\r
 191         }\r
 192 \r
 193         for(int n = 0; n < dimen ; n++){\r
 194             nos[n] = n;\r
 195         }\r
 196 \r
 197         /* nyatla memo\r
 198          * ipが定まらないで計算が行われる場合があるので挿入。\r
 199          * ループ内で0初期化していいかが判らない。\r
 200          */\r
 201         ip=0;\r
 202         //For順変更禁止\r
 203         for(int n=0; n<dimen;n++)\r
 204         {\r
 205             ap_n =ap[n];//wcp = ap + n * rowa;\r
 206             p=0.0;\r
 207             for(int i = n; i<dimen ; i++){//for(i = n, wap = wcp, p = 0.0; i < dimen ; i++, wap += rowa)\r
 208                 if( p < ( pbuf = Math.abs(ap[i][0]))) {\r
 209                     p = pbuf;\r
 210                     ip = i;\r
 211                 }\r
 212             }\r
 213 //          if (p <= matrixSelfInv_epsl){\r
 214             if(p==0.0){\r
 215                 return false;\r
 216 //                throw new NyARException();\r
 217             }\r
 218 \r
 219             nwork  = nos[ip];\r
 220             nos[ip]= nos[n];\r
 221             nos[n] = nwork;\r
 222             \r
 223             ap_ip=ap[ip];\r
 224             for(j=0; j< dimen ; j++){//for(j = 0, wap = ap + ip * rowa, wbp = wcp; j < dimen ; j++) {\r
 225                 work = ap_ip[j];               //work = *wap;\r
 226                 ap_ip[j]=ap_n[j];\r
 227                 ap_n[j]=work;\r
 228             }\r
 229             \r
 230             work=ap_n[0];\r
 231             for(j = 0; j < dimen_1 ; j++){//for(j = 1, wap = wcp, work = *wcp; j < dimen ; j++, wap++)\r
 232                 ap_n[j]=ap_n[j+1]/work;//*wap = *(wap + 1) / work;\r
 233             }\r
 234             ap_n[j]=1.0/work;//*wap = 1.0 / work;\r
 235             for(int i = 0; i < dimen ; i++) {\r
 236                 if(i != n) {\r
 237                     ap_i =ap[i];//wap = ap + i * rowa;\r
 238 \r
 239                     work=ap_i[0];\r
 240                     for(j = 0;j < dimen_1 ; j++){//for(j = 1, wbp = wcp, work = *wap;j < dimen ; j++, wap++, wbp++)\r
 241                         ap_i[j]=ap_i[j+1]-work*ap_n[j];//wap = *(wap + 1) - work * (*wbp);\r
 242                     }\r
 243                     ap_i[j]=-work*ap_n[j];//*wap = -work * (*wbp);\r
 244                 }\r
 245             }\r
 246         }\r
 247 \r
 248         for(int n = 0; n < dimen ; n++) {\r
 249             for(j = n; j < dimen ; j++){\r
 250                 if( nos[j] == n){\r
 251                     break;\r
 252                 }\r
 253             }\r
 254             nos[j] = nos[n];\r
 255             for(int i = 0; i < dimen ;i++){//for(i = 0, wap = ap + j, wbp = ap + n; i < dimen ;i++, wap += rowa, wbp += rowa) {\r
 256                 ap_i=ap[i];\r
 257                 work  =ap_i[j];//work = *wap;\r
 258                 ap_i[j]=ap_i[n];//*wap = *wbp;\r
 259                 ap_i[n]=work;//*wbp = work;\r
 260             }\r
 261         }\r
 262         return true;\r
 263     }\r
 264     /**\r
 265      * sourceの転置行列をdestに得る。arMatrixTrans()の代替品\r
 266      * @param dest\r
 267      * @param source\r
 268      * @return\r
 269      */\r
 270     public static void matrixTrans(NyARMat dest,NyARMat source) throws NyARException\r
 271     {\r
 272         if(dest.row != source.clm || dest.clm != source.row){\r
 273             throw new NyARException();\r
 274         }\r
 275         NyARException.trap("未チェックのパス");\r
 276         //For順変更禁止\r
 277         for(int r=0;r< dest.row;r++){\r
 278             for(int c=0;c<dest.clm;c++){\r
 279                 dest.m[r][c]=source.m[c][r];\r
 280             }\r
 281         }\r
 282     }\r
 283     /**\r
 284      * unitを単位行列に初期化する。arMatrixUnitの代替品\r
 285      * @param unit\r
 286      */\r
 287     public static void matrixUnit(NyARMat unit) throws NyARException\r
 288     {\r
 289         if(unit.row != unit.clm){\r
 290             throw new NyARException();\r
 291         }\r
 292         NyARException.trap("未チェックのパス");\r
 293         //For順変更禁止\r
 294         for(int r = 0; r < unit.getRow(); r++) {\r
 295             for(int c = 0; c < unit.getClm(); c++) {\r
 296                 if(r == c) {\r
 297                     unit.m[r][c]=1.0;\r
 298                 }else{\r
 299                     unit.m[r][c]=0.0;\r
 300                 }\r
 301             }\r
 302         }\r
 303     }\r
 304     /**\r
 305      * sourceの内容を自身に複製する。\r
 306      * Optimized 2008.04.19\r
 307      * @param i_source\r
 308      * @return\r
 309      */\r
 310     public void matrixDup(NyARMat i_source) throws NyARException\r
 311     {\r
 312         //自身の配列サイズを相手のそれより大きいことを保障する。\r
 313         this.realloc(i_source.row,i_source.clm);\r
 314         //内容を転写\r
 315         int r,c;\r
 316         double[][] src_m,dest_m;\r
 317         src_m=i_source.m;\r
 318         dest_m=this.m;\r
 319         //コピーはFor順を変えてもOK\r
 320         for(r = this.row-1; r>=0; r--){\r
 321             for(c =this.clm-1;c>=0; c--)\r
 322             {\r
 323                 dest_m[r][c]=src_m[r][c];\r
 324             }\r
 325         }\r
 326     }\r
 327     public NyARMat matrixAllocDup() throws NyARException\r
 328     {\r
 329         NyARMat result=new NyARMat(this.row,this.clm);\r
 330         //コピー\r
 331         int r,c;\r
 332         double[][] dest_m,src_m;\r
 333         dest_m=result.m;\r
 334         src_m =this.m;\r
 335         //コピーはFor順を変えてもOK\r
 336         for(r = this.row-1; r>=0; r--){\r
 337             for(c =this.clm-1;c>=0; c--)\r
 338             {\r
 339                 dest_m[r][c]=src_m[r][c];\r
 340             }\r
 341         }\r
 342         return result;\r
 343     }    \r
 344     /**\r
 345      * arMatrixInv関数の代替品です。\r
 346      * destにsourceの逆行列を返します。\r
 347      * @param dest\r
 348      * @param source\r
 349      * @throws NyARException\r
 350      */\r
 351     public static void matrixInv(NyARMat dest,NyARMat source) throws NyARException\r
 352     {\r
 353         NyARException.trap("未チェックのパス");\r
 354         dest.matrixDup(source);\r
 355 \r
 356         NyARException.trap("未チェックのパス");\r
 357         dest.matrixSelfInv();\r
 358     }\r
 359     public NyARMat matrixAllocInv() throws NyARException\r
 360     {\r
 361         NyARException.trap("未チェックのパス");\r
 362         NyARMat result=matrixAllocDup();\r
 363 \r
 364         NyARException.trap("未チェックのパス");\r
 365         result.matrixSelfInv();\r
 366         return result;\r
 367     }\r
 368     /**\r
 369      * dim x dim の単位行列を作る。\r
 370      * @param dim\r
 371      * @return\r
 372      * @throws NyARException\r
 373      */\r
 374     public static NyARMat matrixAllocUnit(int dim) throws NyARException\r
 375     {\r
 376         NyARException.trap("未チェックのパス");\r
 377         NyARMat result = new NyARMat(dim, dim);\r
 378         NyARException.trap("未チェックのパス");\r
 379         NyARMat.matrixUnit(result);\r
 380         return result;\r
 381     }\r
 382     /**\r
 383      * arMatrixDispの代替品\r
 384      * @param m\r
 385      * @return\r
 386      */\r
 387     public int matrixDisp() throws NyARException\r
 388     {\r
 389         NyARException.trap("未チェックのパス");\r
 390         System.out.println(" === matrix ("+row+","+clm+") ===");//printf(" === matrix (%d,%d) ===\n", m->row, m->clm);\r
 391         for(int r = 0; r < row; r++){//for(int r = 0; r < m->row; r++) {\r
 392         System.out.print(" |");//printf(" |");\r
 393         for(int c = 0; c < clm; c++) {//for(int c = 0; c < m->clm; c++) {\r
 394             System.out.print(" "+m[r][c]);//printf(" %10g", ARELEM0(m, r, c));\r
 395         }\r
 396         System.out.println(" |");//printf(" |\n");\r
 397         }\r
 398         System.out.println(" ======================");//printf(" ======================\n");\r
 399         return 0;\r
 400     }\r
 401     private final static double PCA_EPS=1e-6;           //#define     EPS             1e-6\r
 402     private final static int            PCA_MAX_ITER=100;       //#define     MAX_ITER        100\r
 403     private final static double PCA_VZERO=1e-16;        //#define     VZERO           1e-16\r
 404     /**\r
 405      * static int EX( ARMat *input, ARVec *mean )の代替関数\r
 406      * Optimize:STEP:[144->110]\r
 407      * @param input\r
 408      * @param mean\r
 409      * @return\r
 410      * @throws NyARException\r
 411      */\r
 412     private void PCA_EX(NyARVec mean) throws NyARException\r
 413     {\r
 414         int lrow,lclm;\r
 415         int i,i2;\r
 416         lrow = this.row;\r
 417         lclm = this.clm;\r
 418         double lm[][]=this.m;\r
 419         \r
 420         if(lrow <= 0 || lclm <= 0){\r
 421             throw new NyARException();\r
 422         }\r
 423         if( mean.getClm() != lclm ){\r
 424             throw new NyARException();\r
 425         }\r
 426 //      double[] mean_array=mean.getArray();\r
 427 //      mean.zeroClear();\r
 428         final double[] mean_array=mean.getArray();\r
 429         double w;\r
 430         //For順変更禁止\r
 431         for(i2=0;i2<lclm;i2++){\r
 432             w=0.0;\r
 433             for(i=0;i<lrow;i++){\r
 434                 //*(v++) += *(m++);\r
 435                 w+=lm[i][i2];\r
 436             }\r
 437             mean_array[i2]=w/lrow;//mean->v[i] /= row;\r
 438         }\r
 439     }\r
 440     /**\r
 441      * static int CENTER( ARMat *inout, ARVec *mean )の代替関数\r
 442      * @param inout\r
 443      * @param mean\r
 444      * @return\r
 445      */\r
 446     private static void PCA_CENTER(NyARMat inout, NyARVec mean) throws NyARException\r
 447     {\r
 448         double[] v;\r
 449         int     row, clm;\r
 450         \r
 451         row = inout.getRow();\r
 452         clm = inout.getClm();\r
 453         if(mean.getClm()!= clm){\r
 454             throw new NyARException();\r
 455         }\r
 456         double[][] im=inout.m;\r
 457         double[] im_i;\r
 458         double w0,w1;\r
 459         v = mean.getArray();\r
 460         //特にパフォーマンスが劣化するclm=1と2ときだけ、別パスで処理します。\r
 461         switch(clm){\r
 462         case 1:\r
 463             w0=v[0];\r
 464             for(int i = 0; i < row; i++ ){\r
 465                 im[i][0]-=w0;\r
 466             }\r
 467             break;\r
 468         case 2:\r
 469             w0=v[0];\r
 470             w1=v[1];\r
 471             for(int i = 0; i < row; i++ ){\r
 472                 im_i=im[i];\r
 473                 im_i[0]-=w0;\r
 474                 im_i[1]-=w1;\r
 475             }\r
 476             break;\r
 477         default:\r
 478             for(int i = 0; i < row; i++ ){\r
 479                 im_i=im[i];\r
 480                 for(int j = 0; j < clm; j++ ){\r
 481                     //*(m++) -= *(v++);\r
 482                     im_i[j]-=v[j];\r
 483                 }\r
 484             }\r
 485         }\r
 486         return;\r
 487     }\r
 488     /**\r
 489      * int x_by_xt( ARMat *input, ARMat *output )の代替関数\r
 490      * @param input\r
 491      * @param output\r
 492      * @throws NyARException\r
 493      */\r
 494     private static void PCA_x_by_xt( NyARMat input, NyARMat output) throws NyARException\r
 495     {\r
 496         NyARException.trap("動作未チェック/配列化未チェック");\r
 497         int     row, clm;\r
 498 //        double[][] out;\r
 499         double[] in1,in2;\r
 500         \r
 501         NyARException.trap("未チェックのパス");\r
 502         row = input.row;\r
 503         clm = input.clm;\r
 504         NyARException.trap("未チェックのパス");\r
 505         if( output.row != row || output.clm != row ){\r
 506             throw new NyARException();\r
 507         }\r
 508         \r
 509 //        out = output.getArray();\r
 510         for(int i = 0; i < row; i++ ) {\r
 511             for(int j = 0; j < row; j++ ) {\r
 512                 if( j < i ) {\r
 513                     NyARException.trap("未チェックのパス");\r
 514                     output.m[i][j]=output.m[j][i];//*out = output->m[j*row+i];\r
 515                 }else{\r
 516                     NyARException.trap("未チェックのパス");\r
 517                     in1=input.m[i];//input.getRowArray(i);//in1 = &(input->m[clm*i]);\r
 518                     in2=input.m[j];//input.getRowArray(j);//in2 = &(input->m[clm*j]);\r
 519                     output.m[i][j]=0;//*out = 0.0;\r
 520                     for(int k = 0; k < clm; k++ ){\r
 521                         output.m[i][j]+=(in1[k]*in2[k]);//*out += *(in1++) * *(in2++);\r
 522                     }\r
 523                 }\r
 524         //                  out.incPtr();\r
 525             }\r
 526         }\r
 527     }\r
 528     /**\r
 529      * static int xt_by_x( ARMat *input, ARMat *output )の代替関数\r
 530      * Optimize:2008.04.19\r
 531      * @param input\r
 532      * @param i_output\r
 533      * @throws NyARException\r
 534      */\r
 535     private static void PCA_xt_by_x(NyARMat input, NyARMat i_output) throws NyARException\r
 536     {\r
 537         double[] in;\r
 538         int     row, clm;\r
 539     \r
 540         row = input.row;\r
 541         clm = input.clm;\r
 542         if(i_output.row!= clm || i_output.clm != clm ){\r
 543             throw new NyARException();\r
 544         }\r
 545         \r
 546         int k,j;\r
 547         double[][] out_m=i_output.m;\r
 548         double w;\r
 549         for(int i = 0; i < clm; i++ ) {\r
 550             for(j = 0; j < clm; j++ ) {\r
 551                 if( j < i ) {\r
 552                     out_m[i][j]=out_m[j][i];//*out = output->m[j*clm+i];\r
 553                 }else{\r
 554                     w=0.0;//*out = 0.0;\r
 555                     for(k = 0; k < row; k++ ){\r
 556                         in=input.m[k];//in=input.getRowArray(k);\r
 557                         w+=(in[i]*in[j]);//*out += *in1 * *in2;\r
 558                     }\r
 559                     out_m[i][j]=w;\r
 560                 }\r
 561             }\r
 562         }\r
 563     }\r
 564     private final NyARVec wk_PCA_QRM_ev=new NyARVec(1);\r
 565     /**\r
 566      * static int QRM( ARMat *a, ARVec *dv )の代替関数\r
 567      * @param a\r
 568      * @param dv\r
 569      * @throws NyARException\r
 570      */\r
 571     private void PCA_QRM(NyARVec dv) throws NyARException\r
 572     {\r
 573         double  w, t, s, x, y, c;\r
 574         int     dim, iter;\r
 575         double[] dv_array=dv.getArray();\r
 576 \r
 577         dim = this.row;\r
 578         if( dim != this.clm || dim < 2 ){\r
 579             throw new NyARException();\r
 580         }\r
 581         if( dv.getClm() != dim ){\r
 582             throw new NyARException();\r
 583         }\r
 584 \r
 585         NyARVec ev = this.wk_PCA_QRM_ev;\r
 586         ev.realloc(dim);\r
 587         double[] ev_array=ev.getArray();\r
 588         if( ev == null ){\r
 589             throw new NyARException();\r
 590         }\r
 591         final double[][] L_m=this.m;\r
 592         this.vecTridiagonalize(dv,ev,1);\r
 593 \r
 594         ev_array[0]=0.0;//ev->v[0] = 0.0;\r
 595         for(int h = dim-1; h > 0; h-- ) {\r
 596             int j = h;\r
 597             while(j>0 && Math.abs(ev_array[j]) > PCA_EPS*(Math.abs(dv_array[j-1])+Math.abs(dv_array[j]))){// while(j>0 && fabs(ev->v[j]) > EPS*(fabs(dv->v[j-1])+fabs(dv->v[j]))) j--;\r
 598                 j--;\r
 599             }\r
 600             if( j == h ){\r
 601                 continue;\r
 602             }\r
 603             iter = 0;\r
 604             do{\r
 605                 iter++;\r
 606                 if( iter > PCA_MAX_ITER ){\r
 607                     break;\r
 608                 }\r
 609                 w = (dv_array[h-1] - dv_array[h]) / 2;//w = (dv->v[h-1] - dv->v[h]) / 2;//ここ？\r
 610                 t = ev_array[h] * ev_array[h];//t = ev->v[h] * ev->v[h];\r
 611                 s = Math.sqrt(w*w+t);\r
 612                 if( w < 0 ){\r
 613                     s = -s;\r
 614                 }\r
 615                 x=dv_array[j] - dv_array[h] + t/(w+s);//x = dv->v[j] - dv->v[h] + t/(w+s);\r
 616                 y=ev_array[j+1];//y = ev->v[j+1];\r
 617                 for(int k = j; k < h; k++ ){\r
 618                     if( Math.abs(x) >= Math.abs(y)){\r
 619                         if( Math.abs(x) > PCA_VZERO ) {\r
 620                             t = -y / x;\r
 621                             c = 1 / Math.sqrt(t*t+1);\r
 622                             s = t * c;\r
 623                         }else{\r
 624                             c = 1.0;\r
 625                             s = 0.0;\r
 626                         }\r
 627                     }else{\r
 628                         t = -x / y;\r
 629                         s = 1.0 / Math.sqrt(t*t+1);\r
 630                         c = t * s;\r
 631                     }\r
 632                     w = dv_array[k] - dv_array[k+1];//w = dv->v[k] - dv->v[k+1];\r
 633                     t = (w * s + 2 * c * ev_array[k+1]) * s;//t = (w * s + 2 * c * ev->v[k+1]) * s;\r
 634                     dv_array[k]-=t;//dv->v[k]   -= t;\r
 635                     dv_array[k+1]+=t;//dv->v[k+1] += t;\r
 636                     if( k > j){\r
 637                         NyARException.trap("未チェックパス");{\r
 638                             ev_array[k]=c * ev_array[k] - s * y;//ev->v[k] = c * ev->v[k] - s * y;\r
 639                         }\r
 640                     }\r
 641                     ev_array[k+1]+=s * (c * w - 2 * s * ev_array[k+1]);//ev->v[k+1] += s * (c * w - 2 * s * ev->v[k+1]);\r
 642 \r
 643                     for(int i = 0; i < dim; i++ ){\r
 644                         x = L_m[k][i];//x = a->m[k*dim+i];\r
 645                         y = L_m[k+1][i];//y = a->m[(k+1)*dim+i];\r
 646                         L_m[k][i]=c * x - s * y;//a->m[k*dim+i]     = c * x - s * y;\r
 647                         L_m[k+1][i]=s * x + c * y;//a->m[(k+1)*dim+i] = s * x + c * y;\r
 648                     }\r
 649                     if( k < h-1 ) {\r
 650                         NyARException.trap("未チェックパス");{\r
 651                             x = ev_array[k+1];//x = ev->v[k+1];\r
 652                             y =-s*ev_array[k+2];//y = -s * ev->v[k+2];\r
 653                             ev_array[k+2]*=c;//ev->v[k+2] *= c;\r
 654                         }\r
 655                     }\r
 656                 }\r
 657             }while(Math.abs(ev_array[h]) > PCA_EPS*(Math.abs(dv_array[h-1])+Math.abs(dv_array[h])));\r
 658         }\r
 659         for(int k = 0; k < dim-1; k++ ) {\r
 660             int h = k;\r
 661             t=dv_array[h];//t = dv->v[h];\r
 662             for(int i = k+1; i < dim; i++ ){\r
 663                 if(dv_array[i] > t ){//if( dv->v[i] > t ) {\r
 664                     h = i;\r
 665                     t=dv_array[h];//t = dv->v[h];\r
 666                 }\r
 667             }\r
 668             dv_array[h]=dv_array[k];//dv->v[h] = dv->v[k];\r
 669             dv_array[k]=t;//dv->v[k] = t;\r
 670             this.flipRow(h,k);\r
 671         }\r
 672     }\r
 673     /**\r
 674      * i_row_1番目の行と、i_row_2番目の行を入れ替える。\r
 675      * @param i_row_1\r
 676      * @param i_row_2\r
 677      */\r
 678     private void flipRow(int i_row_1,int i_row_2)\r
 679     {\r
 680         int i;\r
 681         double w;\r
 682         double[] r1=this.m[i_row_1],r2=this.m[i_row_2];\r
 683         //For順変更OK\r
 684         for(i=clm-1;i>=0;i--){\r
 685             w=r1[i];\r
 686             r1[i]=r2[i];\r
 687             r2[i]=w;\r
 688         }\r
 689     }\r
 690     /**\r
 691      * static int EV_create( ARMat *input, ARMat *u, ARMat *output, ARVec *ev )の代替関数\r
 692      * @param input\r
 693      * @param u\r
 694      * @param output\r
 695      * @param ev\r
 696      * @throws NyARException\r
 697      */\r
 698     private static void PCA_EV_create(NyARMat input, NyARMat u, NyARMat output, NyARVec ev) throws NyARException\r
 699     {\r
 700         NyARException.trap("未チェックのパス");\r
 701         int     row, clm;\r
 702         row = input.row;//row = input->row;\r
 703         clm = input.clm;//clm = input->clm;\r
 704         if( row <= 0 || clm <= 0 ){\r
 705             throw new NyARException();\r
 706        }\r
 707         if( u.row != row || u.clm != row ){//if( u->row != row || u->clm != row ){\r
 708             throw new NyARException();\r
 709         }\r
 710         if( output.row != row || output.clm != clm ){//if( output->row != row || output->clm != clm ){\r
 711             throw new NyARException();\r
 712         }\r
 713         if( ev.getClm()!= row ){//if( ev->clm != row ){\r
 714             throw new NyARException();\r
 715         }\r
 716         double[][] m,in;\r
 717         double[]  m1,ev_array;\r
 718         double  sum, work;\r
 719     \r
 720         NyARException.trap("未チェックのパス");\r
 721         m =output.m;//m = output->m;\r
 722         in=input.m;\r
 723         int i;\r
 724         ev_array=ev.getArray();\r
 725         for(i = 0; i < row; i++ ) {\r
 726             NyARException.trap("未チェックのパス");\r
 727             if( ev_array[i]<PCA_VZERO ){//if( ev->v[i] < VZERO ){\r
 728                 break;\r
 729             }\r
 730             NyARException.trap("未チェックのパス");\r
 731             work = 1 / Math.sqrt(Math.abs(ev_array[i]));//work = 1 / sqrt(fabs(ev->v[i]));\r
 732             for(int j = 0; j < clm; j++ ) {\r
 733                 sum = 0.0;\r
 734                 m1=u.m[i];//m1 = &(u->m[i*row]);\r
 735     //              m2=input.getPointer(j);//m2 = &(input->m[j]);\r
 736                 for(int k = 0; k < row; k++ ) {\r
 737                     sum+=m1[k]+in[k][j];//sum += *m1 * *m2;\r
 738     //                  m1.incPtr();   //m1++;\r
 739     //                  m2.addPtr(clm);//m2 += clm;\r
 740                 }\r
 741                 m1[j]=sum * work;//*(m++) = sum * work;\r
 742     //          {//*(m++) = sum * work;\r
 743     //          m.set(sum * work);\r
 744     //          m.incPtr();}\r
 745             }\r
 746         }\r
 747         for( ; i < row; i++ ) {\r
 748         NyARException.trap("未チェックのパス");\r
 749             ev_array[i]=0.0;//ev->v[i] = 0.0;\r
 750             for(int j = 0; j < clm; j++ ){\r
 751                 m[i][j]=0.0;\r
 752     //          m.set(0.0);//*(m++) = 0.0;\r
 753     //          m.incPtr();\r
 754             }\r
 755         }\r
 756     }\r
 757     private NyARMat wk_PCA_PCA_u=null;\r
 758     /**\r
 759      * static int PCA( ARMat *input, ARMat *output, ARVec *ev )\r
 760      * \r
 761      * @param output\r
 762      * @param o_ev\r
 763      * @throws NyARException\r
 764      */\r
 765     private void PCA_PCA(NyARMat o_output, NyARVec o_ev) throws NyARException\r
 766     {\r
 767         \r
 768         int l_row, l_clm, min;\r
 769         double[] ev_array=o_ev.getArray();\r
 770 \r
 771         l_row =this.row;//row = input->row;\r
 772         l_clm =this.clm;//clm = input->clm;\r
 773         min =(l_clm < l_row)? l_clm: l_row;\r
 774         if( l_row < 2 || l_clm < 2 ){\r
 775             throw new NyARException();\r
 776         }\r
 777         if( o_output.clm != this.clm){//if( output->clm != input->clm ){\r
 778             throw new NyARException();\r
 779         }\r
 780         if( o_output.row!= min ){//if( output->row != min ){\r
 781             throw new NyARException();\r
 782         }\r
 783         if( o_ev.getClm() != min ){//if( ev->clm != min ){\r
 784             throw new NyARException();\r
 785         }\r
 786         \r
 787         NyARMat u;//        u =new NyARMat( min, min );\r
 788         if(this.wk_PCA_PCA_u==null){\r
 789             u=new NyARMat( min, min );\r
 790             this.wk_PCA_PCA_u=u;\r
 791         }else{\r
 792             u=this.wk_PCA_PCA_u;\r
 793             u.realloc(min,min);\r
 794         }\r
 795         \r
 796     \r
 797         if( l_row < l_clm ){\r
 798             NyARException.trap("未チェックのパス");\r
 799             PCA_x_by_xt( this, u );//if(x_by_xt( input, u ) < 0 ) {\r
 800         }else{\r
 801             PCA_xt_by_x( this, u );//if(xt_by_x( input, u ) < 0 ) {\r
 802         }\r
 803         u.PCA_QRM(o_ev);\r
 804 \r
 805         double[][] m1,m2;\r
 806         if( l_row < l_clm ) {\r
 807             NyARException.trap("未チェックのパス");\r
 808             PCA_EV_create( this, u, o_output, o_ev );\r
 809         }else{\r
 810             m1=u.m;//m1 = u->m;\r
 811             m2=o_output.m;//m2 = output->m;\r
 812             int i;\r
 813             for(i = 0; i < min; i++){\r
 814                 if( ev_array[i] < PCA_VZERO){//if( ev->v[i] < VZERO ){\r
 815                     break;\r
 816                 }\r
 817                 for(int j = 0; j < min; j++ ){\r
 818                     m2[i][j]=m1[i][j];//*(m2++) = *(m1++);\r
 819                 }\r
 820             }\r
 821             for( ; i < min; i++){//for( ; i < min; i++){\r
 822                 //コードを見た限りあってそうだからコメントアウト(2008/03/26)NyARException.trap("未チェックのパス");\r
 823                 ev_array[i]=0.0;//ev->v[i] = 0.0;\r
 824                 for(int j = 0; j < min; j++ ){\r
 825                     m2[i][j]=0.0;//*(m2++) = 0.0;\r
 826                 }\r
 827             }\r
 828         }\r
 829     }\r
 830     private NyARMat wk_work_matrixPCA=null;\r
 831     /**\r
 832      * int    arMatrixPCA( ARMat *input, ARMat *evec, ARVec *ev, ARVec *mean );\r
 833      * 関数の置き換え。input引数がthisになる。\r
 834      * Optimize:2008.04.19\r
 835      * @param o_evec\r
 836      * @param o_ev\r
 837      * \r
 838      * @param mean\r
 839      * @throws NyARException\r
 840      */    \r
 841     public void matrixPCA(NyARMat o_evec, NyARVec o_ev, NyARVec mean) throws NyARException\r
 842     {\r
 843         double srow, sum;\r
 844         int     l_row, l_clm;\r
 845         int     check;\r
 846         \r
 847 \r
 848         l_row=this.row;//row = input->row;\r
 849         l_clm=this.clm;//clm = input->clm;\r
 850         check = (l_row < l_clm)? l_row: l_clm;\r
 851         if( l_row < 2 || l_clm < 2 ){\r
 852             throw new NyARException();\r
 853         }\r
 854         if( o_evec.clm != l_clm || o_evec.row != check ){//if( evec->clm != input->clm || evec->row != check ){\r
 855             throw new NyARException();\r
 856         }\r
 857         if( o_ev.getClm()   != check ){//if( ev->clm   != check ){\r
 858             throw new NyARException();\r
 859         }\r
 860         if( mean.getClm() != l_clm){//if( mean->clm != input->clm ){\r
 861             throw new NyARException();\r
 862         }\r
 863 \r
 864         //自分の内容をワークにコピー(高速化の為に、1度作ったインスタンスは使いまわす)\r
 865         NyARMat work;\r
 866         if(this.wk_work_matrixPCA==null){\r
 867             work=this.matrixAllocDup();\r
 868             this.wk_work_matrixPCA=work;\r
 869         }else{\r
 870             work=this.wk_work_matrixPCA;\r
 871             work.matrixDup(this);//arMatrixAllocDup( input );work = arMatrixAllocDup( input );\r
 872         }\r
 873         \r
 874     \r
 875         srow = Math.sqrt((double)l_row);\r
 876         work.PCA_EX(mean );\r
 877 \r
 878         PCA_CENTER(work,mean);\r
 879 \r
 880 \r
 881         int i,j;\r
 882         //For順変更OK\r
 883         for(i=0; i<l_row; i++){\r
 884             for(j=0;j<l_clm;j++){\r
 885                 work.m[i][j]/=srow;//work->m[i] /= srow;\r
 886             }\r
 887         }\r
 888     \r
 889         work.PCA_PCA(o_evec, o_ev);\r
 890     \r
 891         sum = 0.0;\r
 892         double[] ev_array=o_ev.getArray();\r
 893         int ev_clm=o_ev.getClm();\r
 894         //For順変更禁止\r
 895         for(i=0;i<ev_clm;i++){//for(int i = 0; i < ev->clm; i++ ){\r
 896                 sum+=ev_array[i];//sum += ev->v[i];\r
 897         }\r
 898         //For順変更禁止\r
 899         for(i=0;i<ev_clm;i++){//for(int i = 0; i < ev->clm; i++ ){\r
 900                 ev_array[i]/=sum;//ev->v[i] /= sum;\r
 901         }\r
 902     }\r
 903 \r
 904     /*int    arMatrixPCA2( ARMat *input, ARMat *evec, ARVec *ev );*/\r
 905     public static void arMatrixPCA2( NyARMat input, NyARMat evec, NyARVec ev) throws NyARException\r
 906     {\r
 907         NyARException.trap("未チェックのパス");\r
 908         NyARMat   work;\r
 909         // double  srow; // unreferenced\r
 910         double  sum;\r
 911         int     row, clm;\r
 912         int     check;\r
 913 \r
 914         row=input.row;//row = input->row;\r
 915         clm=input.clm;//clm = input->clm;\r
 916         check = (row < clm)? row: clm;\r
 917         if( row < 2 || clm < 2 ){\r
 918             throw new NyARException();\r
 919         }\r
 920         if( evec.getClm()!= input.clm|| evec.row!=check){//if( evec->clm != input->clm || evec->row != check ){\r
 921             throw new NyARException();\r
 922         }\r
 923         if( ev.getClm() != check ){//if( ev->clm   != check ){\r
 924             throw new NyARException();\r
 925         }\r
 926         \r
 927         NyARException.trap("未チェックのパス");\r
 928         work =input.matrixAllocDup();\r
 929 \r
 930         NyARException.trap("未チェックパス");\r
 931         work.PCA_PCA(evec, ev );//rval = PCA( work, evec, ev );\r
 932         sum = 0.0;\r
 933         double[] ev_array=ev.getArray();\r
 934         for(int i = 0; i < ev.getClm(); i++ ){//for( i = 0; i < ev->clm; i++ ){\r
 935             NyARException.trap("未チェックパス");\r
 936             sum+=ev_array[i];//sum += ev->v[i];\r
 937         }\r
 938         for(int i = 0; i < ev.getClm(); i++ ){//for(int i = 0; i < ev->clm; i++ ){\r
 939             NyARException.trap("未チェックパス");\r
 940                 ev_array[i]/=sum;//ev->v[i] /= sum;\r
 941         }\r
 942         return;\r
 943     }\r
 944     public static NyARMat matrixAllocMul(NyARMat a, NyARMat b) throws NyARException\r
 945     {\r
 946         NyARException.trap("未チェックのパス");\r
 947         NyARMat dest=new NyARMat(a.row, b.clm);\r
 948         NyARException.trap("未チェックのパス");\r
 949         dest.matrixMul(a, b);\r
 950         return dest;\r
 951     }\r
 952     /*static double mdet(double *ap, int dimen, int rowa)*/\r
 953     private static double Det_mdet(double[][] ap, int dimen, int rowa) throws NyARException\r
 954     {\r
 955         NyARException.trap("動作未チェック/配列化未チェック");\r
 956         double det = 1.0;\r
 957         double work;\r
 958         int    is = 0;\r
 959         int    mmax;\r
 960     \r
 961         for(int k = 0; k < dimen - 1; k++) {\r
 962             mmax = k;\r
 963             for(int i = k + 1; i < dimen; i++){\r
 964 //              if (Math.abs(arMatrixDet_MATRIX_get(ap, i, k, rowa)) > Math.abs(arMatrixDet_MATRIX_get(ap, mmax, k, rowa))){\r
 965                 if (Math.abs(ap[i][k]) > Math.abs(ap[mmax][k])){\r
 966                     mmax = i;\r
 967                 }\r
 968             }\r
 969             if(mmax != k) {\r
 970                 for (int j = k; j < dimen; j++) {\r
 971                     work = ap[k][j];//work = MATRIX(ap, k, j, rowa);\r
 972                     ap[k][j]=ap[mmax][j];//MATRIX(ap, k, j, rowa) = MATRIX(ap, mmax, j, rowa);\r
 973                     ap[mmax][j]=work;//MATRIX(ap, mmax, j, rowa) = work;\r
 974                 }\r
 975                 is++;\r
 976             }\r
 977             for(int i = k + 1; i < dimen; i++) {\r
 978                 work = ap[i][k]/ ap[k][k];//work = arMatrixDet_MATRIX_get(ap, i, k, rowa) / arMatrixDet_MATRIX_get(ap, k, k, rowa);\r
 979                 for (int j = k + 1; j < dimen; j++){\r
 980                         //MATRIX(ap, i, j, rowa) -= work * MATRIX(ap, k, j, rowa);\r
 981                         ap[i][j]-=work * ap[k][j];\r
 982                 }\r
 983             }\r
 984         }\r
 985         for(int i = 0; i < dimen; i++){\r
 986             det=ap[i][i];//det *= MATRIX(ap, i, i, rowa);\r
 987         }\r
 988         for(int i = 0; i < is; i++){ \r
 989             det *= -1.0;\r
 990         }\r
 991         return det;\r
 992     }\r
 993     /*double arMatrixDet(ARMat *m);*/\r
 994     public static double arMatrixDet(NyARMat m) throws NyARException\r
 995     {\r
 996         NyARException.trap("動作未チェック/配列化未チェック");\r
 997         if(m.row != m.clm){\r
 998             return 0.0;\r
 999         }\r
1000         return Det_mdet(m.getArray(), m.row, m.clm);//return mdet(m->m, m->row, m->row);\r
1001     }\r
1002     private final NyARVec wk_vecTridiagonalize_vec=new NyARVec(0);\r
1003     private final NyARVec wk_vecTridiagonalize_vec2=new NyARVec(0);\r
1004     /**\r
1005      * arVecTridiagonalize関数の代替品\r
1006      * a,d,e間で演算をしてる。何をどうしているかはさっぱりさっぱり\r
1007      * @param a\r
1008      * @param d\r
1009      * @param e\r
1010      * @param i_e_start\r
1011      * 演算開始列(よくわからないけどarVecTridiagonalizeの呼び出し元でなんかしてる)\r
1012      * @return\r
1013      * @throws NyARException\r
1014      */\r
1015     private void vecTridiagonalize(NyARVec d, NyARVec e,int i_e_start) throws NyARException\r
1016     {\r
1017         NyARVec vec=wk_vecTridiagonalize_vec;\r
1018         //double[][] a_array=a.getArray();\r
1019         double  s, t, p, q;\r
1020         int     dim;\r
1021 \r
1022         if(this.clm!=this.row){//if(a.getClm()!=a.getRow()){\r
1023             throw new NyARException();\r
1024         }\r
1025         if(this.clm != d.getClm()){//if(a.getClm() != d.clm){\r
1026             throw new NyARException();\r
1027         }\r
1028         if(this.clm != e.getClm()){//if(a.getClm() != e.clm){\r
1029             throw new NyARException();\r
1030         }\r
1031         dim = this.getClm();\r
1032         \r
1033         double[] d_vec,e_vec;\r
1034         d_vec=d.getArray();\r
1035         e_vec=e.getArray();\r
1036         double[] a_vec_k;\r
1037 \r
1038         for(int k = 0; k < dim-2; k++ ){\r
1039             \r
1040             a_vec_k=this.m[k];\r
1041             vec.setNewArray(a_vec_k,clm);//vec=this.getRowVec(k);//double[] vec_array=vec.getArray();\r
1042             NyARException.trap("未チェックパス"); \r
1043             d_vec[k]=a_vec_k[k];//d.v[k]=vec.v[k];//d.set(k,v.get(k));    //d->v[k] = v[k];\r
1044 \r
1045             //wv1.clm = dim-k-1;\r
1046             //wv1.v = &(v[k+1]);\r
1047             NyARException.trap("未チェックパス"); \r
1048             e_vec[k+i_e_start]=vec.vecHousehold(k+1);//e.v[k+i_e_start]=vec.vecHousehold(k+1);//e->v[k] = arVecHousehold(&wv1);\r
1049             if(e_vec[k+i_e_start]== 0.0 ){//if(e.v[k+i_e_start]== 0.0 ){//if(e.v[k+i_e_start]== 0.0 ){\r
1050                 continue;\r
1051             }\r
1052 \r
1053             for(int i = k+1; i < dim; i++ ){\r
1054                 s = 0.0;\r
1055                 for(int j = k+1; j < i; j++ ) {\r
1056                     NyARException.trap("未チェックのパス");\r
1057                     s += this.m[j][i] * a_vec_k[j];//s += a_array[j][i] * vec.v[j];//s += a.get(j*dim+i) * v.get(j);//s += a->m[j*dim+i] * v[j];\r
1058                 }\r
1059                 for(int j = i; j < dim; j++ ) {\r
1060                     NyARException.trap("未チェックのパス");\r
1061                     s += this.m[i][j] * a_vec_k[j];//s += a_array[i][j] * vec.v[j];//s += a.get(i*dim+j) * v.get(j);//s += a->m[i*dim+j] * v[j];\r
1062                 }\r
1063                 NyARException.trap("未チェックのパス");\r
1064                 d_vec[i]=s;//d.v[i]=s;//d->v[i] = s;\r
1065             }\r
1066         \r
1067 \r
1068             //wv1.clm = wv2.clm = dim-k-1;\r
1069             //wv1.v = &(v[k+1]);\r
1070             //wv2.v = &(d->v[k+1]);\r
1071             a_vec_k=this.m[k];\r
1072             vec.setNewArray(a_vec_k,clm);//vec=this.getRowVec(k);\r
1073 //          vec_array=vec.getArray();\r
1074             NyARException.trap("未チェックパス"); \r
1075             t = vec.vecInnerproduct(d,k+1)/ 2;\r
1076             for(int i = dim-1; i > k; i-- ) {\r
1077                 NyARException.trap("未チェックパス"); \r
1078                 p = a_vec_k[i];//p = v.get(i);//p = v[i];\r
1079                 d_vec[i]-=t*p;q=d_vec[i];//d.v[i]-=t*p;q=d.v[i];//q = d->v[i] -= t*p；\r
1080                 for(int j = i; j < dim; j++ ){\r
1081                     NyARException.trap("未チェックパス"); \r
1082                     this.m[i][j]-=p*(d_vec[j] + q*a_vec_k[j]);//a.m[i][j]-=p*(d.v[j] + q*vec.v[j]);//a->m[i*dim+j] -= p*(d->v[j]) + q*v[j];\r
1083                 }\r
1084             }\r
1085         }\r
1086 \r
1087         if( dim >= 2) {\r
1088             d_vec[dim-2]=this.m[dim-2][dim-2];//d.v[dim-2]=a.m[dim-2][dim-2];//d->v[dim-2] = a->m[(dim-2)*dim+(dim-2)];\r
1089             e_vec[dim-2+i_e_start]=this.m[dim-2][dim-1];//e.v[dim-2+i_e_start]=a.m[dim-2][dim-1];//e->v[dim-2] = a->m[(dim-2)*dim+(dim-1)];\r
1090         }\r
1091 \r
1092         if( dim >= 1 ){\r
1093             d_vec[dim-1]=this.m[dim-1][dim-1];//d.v[dim-1]=a_array[dim-1][dim-1];//d->v[dim-1] = a->m[(dim-1)*dim+(dim-1)];\r
1094         }\r
1095         NyARVec vec2=this.wk_vecTridiagonalize_vec2;\r
1096         for(int k = dim-1; k >= 0; k--) {\r
1097             a_vec_k=this.m[k];\r
1098             vec.setNewArray(a_vec_k,clm);//vec=this.getRowVec(k);//v = a.getPointer(k*dim);//v = &(a->m[k*dim]);\r
1099             if( k < dim-2 ) {\r
1100                 for(int i = k+1; i < dim; i++ ){\r
1101                     //wv1.clm = wv2.clm = dim-k-1;\r
1102                     //wv1.v = &(v[k+1]);\r
1103                     //wv2.v = &(a->m[i*dim+k+1]);\r
1104                     vec2.setNewArray(this.m[i],clm);//vec2=this.getRowVec(i);\r
1105                     \r
1106                     t = vec.vecInnerproduct(vec2,k+1);\r
1107                     for(int j = k+1; j < dim; j++ ){\r
1108                         NyARException.trap("未チェックパス"); \r
1109                         this.m[i][j]-=t*a_vec_k[j];//a_array[i][j]-=t*vec.v[j];//a.subValue(i*dim+j,t*v.get(j));//a->m[i*dim+j] -= t * v[j];\r
1110                     }\r
1111                 }\r
1112             }\r
1113             for(int i = 0; i < dim; i++ ){\r
1114                 a_vec_k[i]=0.0;//v.set(i,0.0);//v[i] = 0.0;\r
1115             }\r
1116             a_vec_k[k]=1;//v.set(k,1);//v[k] = 1;\r
1117         }\r
1118         return;\r
1119     }   \r
1120 }