public static void klu_lsolve(int n, int[] Lip, int Lip_offset,
			int[] Llen, int Llen_offset, double[] LU, int nrhs,
			double[] X, int X_offset)
	{
		double[] x = new double[4] ;
		double lik ;
		/*int[]*/double[] Li ;
		double[] Lx ;
		int k, p, i ;
		int[] len = new int[1] ;
		int[] Li_offset = new int[1] ;
		int[] Lx_offset = new int[1] ;

		switch (nrhs)
		{

			case 1:
				for (k = 0 ; k < n ; k++)
				{
					x [0] = X [X_offset + k] ;
					Li = Lx = GET_POINTER (LU, Lip, Lip_offset, Llen, Llen_offset, Li_offset, Lx_offset, k, len) ;
					/* unit diagonal of L is not stored*/
					for (p = 0 ; p < len[0] ; p++)
					{
						//MULT_SUB (X [Li [p]], Lx [p], x [0]) ;
						X [X_offset + (int) Li [Li_offset[0] + p]] -= Lx [Lx_offset[0] + p] * x [0] ;
					}
				}
				break ;

			case 2:

				for (k = 0 ; k < n ; k++)
				{
					x [0] = X [X_offset + 2*k    ] ;
					x [1] = X [X_offset + 2*k + 1] ;
					Li = Lx = GET_POINTER (LU, Lip, Lip_offset, Llen, Llen_offset, Li_offset, Lx_offset, k, len) ;
					for (p = 0 ; p < len[0] ; p++)
					{
						i = (int) Li [Li_offset[0] + p] ;
						lik = Lx [Lx_offset[0] + p] ;
						//MULT_SUB (X [2*i], lik, x [0]) ;
						X [X_offset + 2*i] -= lik * x [0] ;
						//MULT_SUB (X [2*i + 1], lik, x [1]) ;
						X [X_offset + 2*i + 1] -= lik * x [1] ;
					}
				}
				break ;

			case 3:

				for (k = 0 ; k < n ; k++)
				{
					x [0] = X [X_offset + 3*k    ] ;
					x [1] = X [X_offset + 3*k + 1] ;
					x [2] = X [X_offset + 3*k + 2] ;
					Li = Lx = GET_POINTER (LU, Lip, Lip_offset, Llen, Llen_offset, Li_offset, Lx_offset, k, len) ;
					for (p = 0 ; p < len[0] ; p++)
					{
						i = (int) Li [Li_offset[0] + p] ;
						lik = Lx [Lx_offset[0] + p] ;
						//MULT_SUB (X [3*i], lik, x [0]) ;
						X [X_offset + 3*i] -= lik * x [0] ;
						//MULT_SUB (X [3*i + 1], lik, x [1]) ;
						X [X_offset + 3*i + 1] -= lik * x [1] ;
						//MULT_SUB (X [3*i + 2], lik, x [2]) ;
						X [X_offset + 3*i + 2] -= lik * x [2] ;
					}
				}
				break ;

			case 4:

				for (k = 0 ; k < n ; k++)
				{
					x [0] = X [X_offset + 4*k    ] ;
					x [1] = X [X_offset + 4*k + 1] ;
					x [2] = X [X_offset + 4*k + 2] ;
					x [3] = X [X_offset + 4*k + 3] ;
					Li = Lx = GET_POINTER (LU, Lip, Lip_offset, Llen, Llen_offset, Li_offset, Lx_offset, k, len) ;
					for (p = 0 ; p < len[0] ; p++)
					{
						i = (int) Li [Li_offset[0] + p] ;
						lik = Lx [Lx_offset[0] + p] ;
						//MULT_SUB (X [4*i], lik, x [0]) ;
						X [X_offset + 4*i] -= lik * x [0] ;
						//MULT_SUB (X [4*i + 1], lik, x [1]) ;
						X [X_offset + 4*i + 1] -= lik * x [1] ;
						//MULT_SUB (X [4*i + 2], lik, x [2]) ;
						X [X_offset + 4*i + 2] -= lik * x [2] ;
						//MULT_SUB (X [4*i + 3], lik, x [3]) ;
						X [X_offset + 4*i + 3] -= lik * x [3] ;
					}
				}
				break ;

		}
	}
public static void klu_lsolve ( int n , int [ ] Lip , int Lip_offset , int [ ] Llen , int Llen_offset , double [ ] LU , int nrhs , double [ ] X , int X_offset ) { double [ ] x = new double [ 4 ] ; double lik ; double [ ] Li ; double [ ] Lx ; int k , p , i ; int [ ] len = new int [ 1 ] ; int [ ] Li_offset = new int [ 1 ] ; int [ ] Lx_offset = new int [ 1 ] ; switch ( nrhs ) { case 1 : for ( k = 0 ; k < n ; k ++ ) { x [ 0 ] = X [ X_offset + k ] ; Li = Lx = GET_POINTER ( LU , Lip , Lip_offset , Llen , Llen_offset , Li_offset , Lx_offset , k , len ) ; for ( p = 0 ; p < len [ 0 ] ; p ++ ) { X [ X_offset + ( int ) Li [ Li_offset [ 0 ] + p ] ] -= Lx [ Lx_offset [ 0 ] + p ] * x [ 0 ] ; } } break ; case 2 : for ( k = 0 ; k < n ; k ++ ) { x [ 0 ] = X [ X_offset + 2 * k ] ; x [ 1 ] = X [ X_offset + 2 * k + 1 ] ; Li = Lx = GET_POINTER ( LU , Lip , Lip_offset , Llen , Llen_offset , Li_offset , Lx_offset , k , len ) ; for ( p = 0 ; p < len [ 0 ] ; p ++ ) { i = ( int ) Li [ Li_offset [ 0 ] + p ] ; lik = Lx [ Lx_offset [ 0 ] + p ] ; X [ X_offset + 2 * i ] -= lik * x [ 0 ] ; X [ X_offset + 2 * i + 1 ] -= lik * x [ 1 ] ; } } break ; case 3 : for ( k = 0 ; k < n ; k ++ ) { x [ 0 ] = X [ X_offset + 3 * k ] ; x [ 1 ] = X [ X_offset + 3 * k + 1 ] ; x [ 2 ] = X [ X_offset + 3 * k + 2 ] ; Li = Lx = GET_POINTER ( LU , Lip , Lip_offset , Llen , Llen_offset , Li_offset , Lx_offset , k , len ) ; for ( p = 0 ; p < len [ 0 ] ; p ++ ) { i = ( int ) Li [ Li_offset [ 0 ] + p ] ; lik = Lx [ Lx_offset [ 0 ] + p ] ; X [ X_offset + 3 * i ] -= lik * x [ 0 ] ; X [ X_offset + 3 * i + 1 ] -= lik * x [ 1 ] ; X [ X_offset + 3 * i + 2 ] -= lik * x [ 2 ] ; } } break ; case 4 : for ( k = 0 ; k < n ; k ++ ) { x [ 0 ] = X [ X_offset + 4 * k ] ; x [ 1 ] = X [ X_offset + 4 * k + 1 ] ; x [ 2 ] = X [ X_offset + 4 * k + 2 ] ; x [ 3 ] = X [ X_offset + 4 * k + 3 ] ; Li = Lx = GET_POINTER ( LU , Lip , Lip_offset , Llen , Llen_offset , Li_offset , Lx_offset , k , len ) ; for ( p = 0 ; p < len [ 0 ] ; p ++ ) { i = ( int ) Li [ Li_offset [ 0 ] + p ] ; lik = Lx [ Lx_offset [ 0 ] + p ] ; X [ X_offset + 4 * i ] -= lik * x [ 0 ] ; X [ X_offset + 4 * i + 1 ] -= lik * x [ 1 ] ; X [ X_offset + 4 * i + 2 ] -= lik * x [ 2 ] ; X [ X_offset + 4 * i + 3 ] -= lik * x [ 3 ] ; } } break ; } }
Solve Lx=b.  Assumes L is unit lower triangular and where the unit diagonal
entry is NOT stored.  Overwrites B  with the solution X.  B is n-by-nrhs
and is stored in ROW form with row dimension nrhs.  nrhs must be in the
range 1 to 4.

@param n
@param Lip
@param Llen
@param LU
@param nrhs
@param X right-hand-side on input, solution to Lx=b on output