
  'Example of matrix layout mapped to an array of 4x4 cells
  '
  '  0 4 8 C
  '  1 5 9 D
  '  2 6 A E
  '  3 7 B F
  '

  % MatrixType double

  sub MatrixMul(MatrixType *A,*B,*C, sys n)
  '========================================
  '
  '
  #if typecodeof A=0x64
    % OneStep 4
    % mtype single
  #endif
  '
  #if typecodeof A=0x68
    % OneStep 8
    % mtype double
  #endif

  sys pa=@A, pb=@B, pc=@C
  sys ColStep=OneStep*n

  mov ecx,pa
  mov edx,pb
  mov eax,pc

  mov esi,n
  (
  call column : dec esi : jg repeat
  )
  exit sub

  column:
  '======

  mov edi,n
  (
  call cell : dec edi : jg repeat
  )
  add edx,ColStep
  sub ecx,ColStep
  ret

  cell: ' row A * column B
  '=======================

  'matrix data is stored ascending vertically then horizontally
  'thus rows are minor, columns are major
  '
  push  ecx
  push  edx
  push  eax
  mov   eax,4
  fldz
  (
  fld   mtype [ecx]
  fmul  mtype [edx]
  faddp st1
  add   ecx,ColStep 'next column of matrix A
  add   edx,OneStep 'next row of matrix B
  dec eax
  jnz repeat
  )
  pop eax
  fstp mtype [eax] 'assign to next row of matrix C
  '
  pop edx
  pop ecx
  add   eax,OneStep 'next cell in column of matrix C (columns then rows)
  add   ecx,OneStep 'next row of matrix A
  ret
  '
  end sub


  sub MatrixInv(MatrixType *M, sys n)
  '==================================

  'Uses the Gauss-Jordan Method.
  'Inverts matrix, m[], "in-place".

  'limit to max 16x16 matrix

  MatrixType test,colmax,factor,temp
  sys        pivot,row,col,maxrow
  sys        rowswitch,switchcount
  sys        i,j
  sys        rs[128] 'for limit array n=64

  '
  'TO RESOLVE MATRIX INDEX
  '=======================
  
  indexbase 1
  #define mat(row,col) m[ (col-1)*n+row ]
  #define rowswitches(row,col) rs[ col*64-64+row ]
  '
  '
  'TRAVERSE COLUMNS
  '================


  for pivot=1 to n 'VERTICAL SWEEP
  '
  '
  'GET THE MAXIMUM ELEMENT INTO EACH PIVOT COLUMN
  '==============================================
  '
  colmax=0
  maxrow=pivot
  '
  for row=pivot to n 'HORIZONTAL SWEEP
    test=abs (mat(row,pivot))
    if test>colmax then
      colmax=test 'maximum value for each column
      maxrow=row  'cell of maximum value
    end if
  next
  '
  'NEED TO SWAP ROWS?
  '
  if maxrow<>pivot then
    '
    '
    'SWAP ROWS TO GET THE MAX VALUE INTO THE PIVOT COLUMN
    '
    for col=1 to n 'minor
      temp=mat(pivot,col)
      mat(pivot,col)=mat(maxrow,col)
      mat(maxrow,col)=temp
    next
    '
    '
    'REMEMBER THE ROW SWAP PAIRS
    '
    switchcount++
    rowswitches(switchcount,1)=pivot
    rowswitches(switchcount,2)=maxrow
  end if
  '
  '
  'DIVIDE THE PIVOT COLUMN MEMBERS BY THE PIVOT DIAGONAL
  '
  'AND REPLACE THE PIVOT DIAGONAL WITH 1
  '
  factor=1/mat(pivot,pivot)
  mat(pivot,pivot)=1
  '
  for col=1 to n 'minor
    #show mat(pivot,col) *= factor
  next


  'APPLY PIVOT FACTORS TO ALL ELEMENTS, EXCEPT FOR PIVOT ROW MEMBERS
  '
  'also setting each pivot row element to 0
  '
  for row=1 to n 'major
    if row<>pivot then
      factor=-mat(row,pivot)
      mat(row,pivot)=0
      for col=1 to n 'minor
        mat(row,col) += mat(pivot,col) * factor
      next
    end if
  next

  next 'MOVE PIVOT TO NEXT COLUMN


  '
  'SWAP COLUMNS ACCORDING TO ORIGINAL ROW POSITIONS
  '
  for rowswitch=switchcount to 1 step-1 'EACH ROW SWAP RECORD
    i=rowswitches(rowswitch,1)
    j=rowswitches(rowswitch,2)
    for row=1 to n 'SWAP COLUMN-PAIR MEMBERS
      temp = mat(row,i)
      mat(row,i)=mat(row,j)
      mat(row,j)=temp
    next 
  next
  '
  '
  end sub


  function ShowMatrix(MatrixType*A,sys n) as string
  '================================================
    string cr=chr(13)+chr(10), tab=chr(9)
    function="MATRIX " n "x" n cr cr
    sys i,j,m
    '
    for i=1 to n
      m=0
      for j=1 to n
        function+=str( A[m+i] ) tab
        m+=n
      next
      function+=cr
    next
  end function

  'TEST
  '====       

  % n 4
  MatrixType A[n*n],B[n*n],C[n*n]


  'reading vertically (minor) then left to right (major)

  A <= 4,0,0,1, 0,4,0,0, 0,0,4,0, 0,0,0,4

  B <= 2,0,0,2, 0,2,0,0, 0,0,2,0, 0,0,0,2


  MatrixMul A,B,C,n

  Print ShowMatrix C,n

  MatrixInv C,n

  Print ShowMatrix C,n

  MatrixInv C,n

  Print ShowMatrix C,n


  