def rref(A) { S=size(A); M=S[0]; N=S[1]; I=0; for (J = 0; J < N; J++ ) { P = -1; for ( K = I; K < M; K++ ) if ( P < 0 && A[K][J] != 0 ) P = K; printf("(P,J)=(%a,%a) Pivot position.\n",P,J); if ( P >= 0 ) { if ( P != I ) { for ( L = J; L < N; L++ ) { T = A[I][L]; A[I][L] = A[P][L]; A[P][L] = T; } printf("P=%a, I=%a. Swap rows.\n%a\n",P,I,A); } for ( L = J+1; L < N; L++ ) A[I][L] = A[I][L]/A[I][J]; A[I][J] = 1; for ( K = 0; K < M; K++ ) { if ( K != I ) { for ( L = J+1; L < N; L++ ) A[K][L] = A[K][L] - A[K][J]*A[I][L]; A[K][J] = 0; } } printf("Eliminate by I=%a-th row.\n%a\n",I,A); I++; } } return I; } A=newmat(2,3,[[0,2,3],[4,5,6]]); A0=A*1$ rref(A); print(A)$ end$