Modify sgesvx(GEPP)

Clapack 是個由 c 寫成的矩陣計算的 library，最原始的 lapack 是由 fortran 寫成，都可以在 NETLIB 找到下載點。

一般以直接法解線性方程組時，最直接的想法是使用高斯消去法。但在當矩陣不平衡時，做列運算容易產生不必要的錯誤，為了減少這種情況的發生，我們會做 pivoting，而最常用的有 partial pivoting 以及 complete pivoting，通常是使用 partial pivoting 以便在速度以及精確度上取得較佳的平衡。以 partial pivoting 做高斯消去法便稱做 GEPP。

此處大約記錄一下在 clapack 中，如何參考它的 subroutine sgesvx(GEPP) 來另外建立另一種高斯消去法(GECP)。

首先至 clapack 下載 source code，接著解壓縮，假設解壓縮後的資料夾名字是 CLAPACK:
$ cd CLAPACK/
$ cp make.inc.example make.inc
$ cp SRC/sgesvx.c SRC/gecp.c
$ cp SRC/sgetf2.c SRC/sgetf22.c
$ vi SRC/gecp.c

(-)int sgesvx_(char *fact, char *trans, integer *n, integer *
(-)    nrhs, real *a, integer *lda, real *af, integer *ldaf, integer *ipiv,;
(+)int gecp_(char *fact, char *trans, integer *n, integer *
(+)    nrhs, real *a, integer *lda, real *af, integer *ldaf, integer *ipiv, integer *ipiv2,
...
(-)     real r__1, r__2;
(+)     real r__1, r__2, tmp;
(+)     int i;
...
           --ipiv;
(+)     --ipiv2;
...
(-)      sgetrf_(n, n, &af[af_offset], ldaf, &ipiv[1], info);
(+) sgetf22_(n, n, &af[af_offset], ldaf, &ipiv[1], &ipiv2[1], info);
...
   sgetrs_(trans, n, nrhs, &af[af_offset], ldaf, &ipiv[1], &x[x_offset], ldx, info);
(+)     for (i=(*n-1) ; i>=0 ; --i){
(+)         tmp = x[x_offset+i];
(+)                 x[x_offset+i] = x[x_offset+ipiv2[i+1]-1];
(+)                 x[x_offset+ipiv2[i+1]-1] = tmp;
(+)     }

存檔。接著
$ vi SRC/sgetf22.c

(-) int sgetf2_(integer *m, integer *n, real *a, integer *lda,
(-)     integer *ipiv, integer *info)
(+) int sgetf22_(integer *m, integer *n, real *a, integer *lda,
(+)     integer *ipiv, integer *ipiv2, integer *info)
...
(+)     integer tmp, tmp2, *tmp3, tmp4, jp2, matmax;
(+)     int i, k;
(+)     tmp4 = 1;
(+)     tmp3 = &tmp4;
...
...
   --ipiv;
(+)     --ipiv2;
   /* Find pivot and test for singularity. */
   i__2 = *m - j + 1;
(+)     matmax = isamax_(&i__2, &a[j + j * a_dim1], &c__1);
   jp = j - 1 + matmax;
(+)     matmax = jp + j*a_dim1;
(+)     for ( k=1 ; k<=i__1-j ; ++k ){ (+)             tmp = isamax_(&i__2, &a[j + (j+k) * a_dim1], &c__1); (+)             jp2 = j - 1 + tmp; (+)                 if ( a[jp2+(j+k)*a_dim1] > a[matmax] ){
(+)                     matmax = jp2+(j+k)*a_dim1;
(+)             }
(+)             tmp2 = j-1+tmp;
(+)     }
(-)     ipiv[j] = jp;
(+)     tmp = matmax%a_dim1;
(+)     tmp = (tmp==0?a_dim1:tmp);
(+)     tmp2 = ((matmax - a_offset)/a_dim1)+1;
(+)     ipiv[j] = tmp;
(+)     ipiv2[j] = tmp2;
   if (a[jp + j * a_dim1] != 0.f) {
   /* Apply the interchange to columns 1:N. */
(-)             if (jp != j) {
(+)                 if (tmp != j){
(-)                         sswap_(n, &a[j + a_dim1], lda, &a[jp + a_dim1], lda);
(+)                         sswap_(n, &a[j + a_dim1], lda, &a[tmp + a_dim1], lda);
                    }
(+)                  if (tmp2 != j){
(+)                         sswap_(n, &a[1 + j*a_dim1], tmp3, &a[1 + tmp2*a_dim1], tmp3);
(+)              }
   ...
         }

存檔。再來在 SRC/Makefile 裡頭加入 gecp.o 和 sgetf22.o，存檔。最後在 SRC 下 make，這需要一點時間，最後會產生 lapack_LINUX.a。最後還需要至 BLAS/SRC/ 下去 make 出 blas_LINUX.a 以及 F2CLIBS/libf2c/ 下 make libf2c.a。到此為止即完成了新 subroutine gecp 了。後續的驗證，需要另外寫一個小小程式來呼叫 sgesvx 以及 gecp 這兩個 function，觀查 GEPP 以及 GECP 是否在 solution, backward error 以及 forward error 上有所差異。