rSLZ is a simple but useful tool that I used to help create examples for Linear Algebra classes that I taught over the years. In short, it generates random integer matrices with integer inverses and some measure of control over the sizes of the entries. For more details, see the docstring.

rSLZ.py

"""This is a small tool that I used myself while teaching Linear Algebra classes over the years. I frequently found myself in the following situation: 
- I want a "random-looking" diagonalizable matrix with eigenvalues of my choice (including multiplicities), because I want to illustrate certain behavior in class or test certain ideas on an exam. 
- But I also want the matrix to have small integer entries, and I possibly want to know eigenvectors in advance. 
Fundamentally, all you need to do is conjugate (= "similarity transformation") a given diagonal matrix by a "small" integral matrix g whose inverse is also integral and "small". 
A source of most, but not all, such g is GL(Z), the group of integer matrices whose determinants are either +1 or -1. The subgroup SL(Z) of those with determinant +1 is essentially the same thing, and the function rSLZ below produces random elements g of SL(Z). A typical call for classroom/exam use might be rSLZ(3,9) or rSLZ(4,5).
Besides the size of the matrix, you also provide a bound B on the entries of the returned matrix. It is not 100% true that the entries are always bounded by B, but it is "mostly" true. 
Anyway, lives are not at risk here: if you don't like your matrix, generate another one. 
Not much control is exerted over the entries of the inverse -- despite the bound on g, the inverse h of g can have unacceptably large entries. You may need to check a few g before finding one with an acceptable inverse.
If you need something better than this, let me know and I'll look into it -- rSLZ is not very sophisticated but it was easy to write and it did its job well.
AUTHOR: Sean Rostami
CONTACT: sean.rostami@gmail.com"""


import math, random # not seeded here!

import numpy


def _rSLZb( d, i, B ):
	"""[PRIVATE] Assuming that one wants L*U to have entries bounded by B, returns the bound that L's (i,-) and U's (-,i) entries should satisfy."""
	return int( math.ceil( math.sqrt( (B-1)/float(i) ) ) ) # does a good job overall but does not truly guarantee desired bound on L*U... investigate later


def _rSLZT( dim, B, L, dtype = numpy.intc ):
	"""[PRIVATE] Randomly generates a unipotent triangular matrix T. Whether T is upper- or lower-triangular is decided by the parameter L (True <=> lower). 
	If upper- U and lower- L are generated by this function, it is "mostly" true that all entries of L*U are bounded in magnitude by B. 
	The function uses random.randint, and it is the caller's responsibility to seed the generator."""
	T = numpy.identity( dim, dtype )
	for r in range( dim ):
		for c in range( r ):
			b = _rSLZb( dim, r, B )
			T[ ( (r,c) if L else (c,r) ) ] = random.randint( -b, b ) # includes both endpoints
	return T


def _rSLZL( dim, B, dtype = numpy.intc ):
	"""[PRIVATE] Randomly generates and returns a unipotent lower-triangular matrix L. 
	If U is returned by _rSLZU then it is "mostly" true that all entries of L*U are bounded in magnitude by B. 
	The function uses random.randint, and it is the caller's responsibility to seed the generator."""
	return _rSLZT( dim, B, True, dtype )


def _rSLZU( dim, B, dtype = numpy.intc ):
	"""[PRIVATE] Randomly generates and returns a unipotent upper-triangular matrix U. 
	If L is returned by _rSLZL then it is "mostly" true that all entries of L*U are bounded in magnitude by B. 
	The function uses random.randint, and it is the caller's responsibility to seed the generator."""
	return _rSLZT( dim, B, False, dtype )


def rGLZP( dim, dtype = numpy.intc ):
	"""Returns a random dim x dim permutation matrix. 
	The function uses random.sample, and it is the caller's responsibility to seed the generator."""
	P = numpy.zeros( (dim,dim), dtype )
	p = random.sample( range( dim ), dim ) # use numpy.random.shuffle or numpy.random.permutation instead?
	for (i,j) in enumerate( p ):
		P[i,j] = 1
	return P


def rSLZ( dim, B, dtype = numpy.intc ):
	"""Randomly generates and returns a dim x dim element of SL(Z) whose entries are "mostly" bounded in magnitude by B. 
	The function uses random.randint, and it is the caller's responsibility to seed the generator. 
	Note that floating-point inaccuracy may cause numpy.linalg.det to return a value slightly different from 1. 
	Similarly, it is a mathematical fact that the inverse of a matrix returned by rSLZ is also integral, but numpy.linalg.inv may return a matrix with many slightly non-integral entries. 
	If you want the true inverse, apply numpy.round_ to that inverse (or use rSLZ.SLZi). 
	Via the third (optional) parameter, you can change the type of integer stored by these matrices, but I can't imagine that you would ever need or want to do so."""
	return numpy.matmul( _rSLZL( dim, B, dtype ), _rSLZU( dim, B, dtype ) )


def SLZi( x ):
	"""Convenience: Rounds a matrix that "should" be integral but isn't due to floating-point inaccuracy. Code:
	return numpy.round_( numpy.linalg.inv( x ) )"""
	return numpy.round_( numpy.linalg.inv( x ) )