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