proc invers (matrix A)
"USAGE:      invers(A); matrix A
RETURN:      matrix, die Inverse von A, falls A invertierbar, sonst die Nullmatrix
KEYWORDS:    Matrix, Inverse, Gauss-Algorithmus
EXAMPLE:     example invers;   zeigt ein Beispiel"
{
  int m=nrows(A);     // Anzahl der Zeilen von A
  int n=ncols(A);     // Anzahl der Spalten von A
  int i,j;            // Laufindizes 
  matrix B[m][2*n];   // Hilfsmatrix, in die A und E eingetragen werden, um A^{-1} zu berechnen
  // Falls die Matrix nicht quadratisch ist, gib die Nullmatrix zurueck, ...
  if (m!=n)
  {
    return(matrix(0));
  }
  // ... sonst ...
  else
  {
    // ... fuege A und die Einheitsmatrix E in einem Schema zusammen ...
    for (i=1;i<=m;i++)
    {
      for (j=1;j<=n;j++)
      {
        B[i,j]=A[i,j];
      }
      B[i,i+n]=1;
    }
    // ... berechne die reduzierte Zeilen-Stufen-Form des Schemas B ...
    B=RZSF(B);
    // ... speichere den ersten Teil des Schemas in A ab.
    A=B[1..m,1..n];
    // Ist der erste Teil des Schemas die Einheitsmatrix, so ist A invertierbar, ...
    if (A==unitmat(n))
    {
      // ... und die Inverse steht im zweiten Teil des Schemas, ...
      A=B[1..m,n+1..2*n];
      return(A);
    }
    // ... sonst ist die Matrix nicht invertierbar.
    else
    {
      return(matrix(0));
    }
  }
}
example
{
  "EXAMPLE:";
  echo=2;
  ring r=0,x,lp;    
  matrix A[4][4]=1,3,-1,4,2,5,-1,3,0,4,-3,1,-3,1,-5,-2;
  print(invers(A));
  print(invers(A)*A);
  matrix B[4][4]=1,-2,2,0,2,-7,10,-1,-2,4,-5,2,3,-5,7,1;
  print(invers(B));
  print(invers(B)*B);
  matrix C[3][3]=1,2,3,4,5,6,7,8,9;
  print(C);
  invers(C);
  matrix D[4][3]=1,2,3,4,5,6,7,8,9,10,11,13;
  print(D);
  invers(D);
}