역행렬 (가우스 소거법)

 

 

 

 

Javascript로 구성한 Gauss Jordan Elimination 역행렬 계산 알고리즘

//
// Gauss Jordan Elimination
//

	// 0. set unit matrix
	minv = [];
	for( r = 0; r < mat.length; r++){
		
		minv[r] = [];
		for( c = 0; c < mat[r].length; c++){		
			minv[r][c] = 0;
		}
		
		minv[r][r] = 1;
	}

// 1. calc inv mat

	// Start solving.
	for( r = 0; r < mat.length ; r++){ 								// by row
		// Zero out all entries in column r after this row.
		// See if this row has a non-zero entry in column r.
		if(Math.abs( mat[r][r] ) < Math.pow(1/10,10) ){
			// Zero value.
			// Try to swap with a later row.
			for( r2 = r + 1; r2 < mat.length; r2++){
				if(Math.abs( mat[r2][r] ) > Math.pow(1/10,10) ){
					// This row will work. Swap them.
					for( c = 0; c < mat[r2].length; c++){
						tmp = mat[r][c];
						mat[r][c] = mat[r2][c];
						mat[r2][c] = tmp;
					}
					// This row will work. Swap them.
					for( c = 0; c < minv[r2].length; c++){
						tmp = minv[r][c];
						minv[r][c] = minv[r2][c];
						minv[r2][c] = tmp;
					}
					
					break;
				}
			}
		}
					
		// If this row has a non-zero entry in column r,
		if(Math.abs( mat[r][r] ) > Math.pow(1/10,10) ){
			
			// 현재 열의 첫번째 값이 1이 되도록 연산
			factor = mat[r][r];
			for( c = 0; c < mat[r].length; c++){  // 시작 위치 확인 필요
				mat[r][c] = mat[r][c] / factor ;
				minv[r][c] = minv[r][c] / factor ;
			}
			
			mat[r][r] = 1;				
			
			// 가로 열 계산
			//for( r2 = r+1; r2 < r+2 ; r2++){  // from Top to End row						
			for( r2 = r+1; r2 < mat.length ; r2++){  // from Top to End row						
				factor = -mat[r2][r] / mat[r][r];
				for( c = 0; c < mat[r2].length; c++){  // 시작 위치 확인 필요
					mat[r2][c] = mat[r2][c] + factor * mat[r][c];
					minv[r2][c] = minv[r2][c] + factor * minv[r][c];
				}	
					
			}
			
		}
	}
	
	// backward solve
	for( r = mat.length - 1; r >= 0 ; r--){
		for( r2 = r - 1; r2 >= 0 ; r2--){
			factor = -mat[r2][r] / mat[r][r];
			for( c = 0; c < mat[r].length; c++){
				mat[r2][c] = mat[r2][c] + factor * mat[r][c];
				minv[r2][c] = minv[r2][c] + factor * minv[r][c];
			}
		}
	}		
	
	return minv;	

 

 

 

 

 

 

Excel VBA로 구성한 Gauss Jordan Elimination 역행렬 계산 알고리즘

 

Private Function Inv_Mat(aMat)
'-----------------------------------------------------------------
' Gauss Jordan Elimination
'-----------------------------------------------------------------

Dim aMinv()
Dim ir, ir2, ic As Integer
Dim dTmp, dFactor As Double

' 0. set unit matrix
    ReDim aMinv(UBound(aMat), UBound(aMat, 2))
    
    For ir = 0 To UBound(aMat)
        
        'aminv[r] = [];
        For ic = 0 To UBound(aMat, 2)
        
            aMinv(ir, ic) = 0
            
        Next ic
        
        aMinv(ir, ir) = 1
        
    Next ir

' 1. calc inv mat

    '/ Start solving.
    For ir = 0 To UBound(aMat)                              '/ by row
        '/ Zero out all entries in column r after this row.
        '/ See if this row has a non-zero entry in column r.
        If Math.Abs(aMat(ir, ir)) < (1 / 10) ^ 10 Then
            '/ Zero value.
            '/ Try to swap with a later row.
            '/ 대각행렬의 값이 0이면, 다음줄의 값과 교체..
            For ir2 = ir + 1 To UBound(aMat)
            
                If Math.Abs(aMat(ir2, ir) > 1 / 10 ^ 10) Then
                    '// This row will work. Swap them.
                    For ic = 0 To UBound(aMat, 2)
                        dTmp = aMat(ir, ic)
                        aMat(ir, ic) = aMat(ir2, ic)
                        aMat(ir2, ic) = dTmp
                    Next ic
                    
                    '// This row will work. Swap them.
                    For ic = 0 To UBound(aMinv, 2)
                        dTmp = aMinv(ir, ic)
                        aMinv(ir, ic) = aMinv(ir2, ic)
                        aMinv(ir2, ic) = dTmp
                    Next ic
                    
                    Exit For
                    
                End If
                
            Next ir2
            
        End If
                    
        '// If this row has a non-zero entry in column r,
        If Math.Abs(aMat(ir, ir)) > 1 / 10 ^ 10 Then
            
            '// 현재 열의 첫번째 값이 1이 되도록 연산
            dFactor = aMat(ir, ir)
            
            For ic = 0 To UBound(aMat, 2)   '// 시작 위치 확인 필요
                aMat(ir, ic) = aMat(ir, ic) / dFactor
                aMinv(ir, ic) = aMinv(ir, ic) / dFactor
            Next ic
            
            aMat(ir, ir) = 1
            
            '// 가로 열 계산
            '//for( r2 = r+1; r2 < r+2 ; r2++){  // from Top to End row
            For ir2 = ir + 1 To UBound(aMat)   '/ from Top to End row
                dFactor = -aMat(ir2, ir) / aMat(ir, ir)
                For ic = 0 To UBound(aMat, 2)     '// 시작 위치 확인 필요
                    aMat(ir2, ic) = aMat(ir2, ic) + dFactor * aMat(ir, ic)
                    aMinv(ir2, ic) = aMinv(ir2, ic) + dFactor * aMinv(ir, ic)
                Next ic
                    
            Next ir2
            
        End If
        
    Next ir
    
    '// backward solve
    For ir = UBound(aMat) To 0 Step -1
    
        For ir2 = ir - 1 To 0 Step -1
        
            dFactor = -aMat(ir2, ir) / aMat(ir, ir)
            
            For ic = 0 To UBound(aMat)
            
                aMat(ir2, ic) = aMat(ir2, ic) + dFactor * aMat(ir, ic)
                aMinv(ir2, ic) = aMinv(ir2, ic) + dFactor * aMinv(ir, ic)
                
            Next ic
            
        Next ir2
        
    Next ir
    
' return result
Inv_Mat = aMinv

End Function