There are several restrictions. The first is that there are only a limited number of rings for which this function is implemented. Second, over
RR or
CC, the matrix
A must be a square non-singular matrix. Third, if
A and
b are mutable matrices over
RR or
CC, they must be dense matrices.
i1 : kk = ZZ/101;
|
i2 : A = matrix"1,2,3,4;1,3,6,10;19,7,11,13" ** kk
o2 = | 1 2 3 4 |
| 1 3 6 10 |
| 19 7 11 13 |
3 4
o2 : Matrix kk <--- kk
|
i3 : b = matrix"1;1;1" ** kk
o3 = | 1 |
| 1 |
| 1 |
3 1
o3 : Matrix kk <--- kk
|
i4 : x = solve(A,b)
o4 = | 2 |
| -1 |
| 34 |
| 0 |
4 1
o4 : Matrix kk <--- kk
|
i5 : A*x-b
o5 = 0
3 1
o5 : Matrix kk <--- kk
|
Over
RR or
CC, the matrix
A must be a non-singular square matrix.
i6 : kk = RR
o6 = RR
o6 : Ring
|
i7 : A = matrix"1,2,3;1,3,6;19,7,11" ** kk
o7 = | 1 2.000000 3.000000 |
| 1 3.000000 6.000000 |
| 19.000000 7.000000 11.000000 |
3 3
o7 : Matrix RR <--- RR
|
i8 : b = matrix"1;1;1" ** kk
o8 = | 1 |
| 1 |
| 1 |
3 1
o8 : Matrix RR <--- RR
|
i9 : x = solve(A,b)
o9 = | -0.148936 |
| 1.148936 |
| -0.382979 |
3 1
o9 : Matrix RR <--- RR
|
i10 : A*x-b
o10 = 0
3 1
o10 : Matrix RR <--- RR
|
For large dense matrices over
RR or
CC, this function calls the lapack routines.
i11 : kk = CC
o11 = CC
o11 : Ring
|
i12 : A = mutableZero(kk,1000,1000, Dense=>true)
o12 = 0
o12 : MutableMatrix
|
i13 : b = mutableZero(kk,1000,2, Dense=>true)
o13 = 0
o13 : MutableMatrix
|
i14 : for i from 1 to 10000 do A_(random 1000, random 1000) = random 1.0 + ii * random 1.0;
|
i15 : for i from 0 to 999 do b_(i,0) = random 1.0 + ii * random 1.0;
|
i16 : for i from 0 to 999 do b_(i,1) = random 1.0 + ii * random 1.0;
|
i17 : time x = solve(A,b);
-- used 2.46815 seconds
|
i18 : C = (matrix A)*(matrix x)-matrix b;
1000 2
o18 : Matrix CC <--- CC
|
i19 : C_(123,0)
o19 = - 3.32165*10^-14 + 7.93809*10^-14ii
o19 : CC
|
i20 : C_(567,1)
o20 = - 5.93525*10^-13 - 4.56302*10^-14ii
o20 : CC
|
This may be used to invert a matrix over
ZZ/p,
RR or
QQ.
i21 : kk = RR
o21 = RR
o21 : Ring
|
i22 : A = random(kk^5, kk^5)
o22 = | 0.220789 -0.979717 -0.586282 0.495034 -0.551757 |
| -0.252843 -0.526080 -0.895040 -0.933781 0.234187 |
| 0.395667 -0.476370 -0.157478 0.556580 0.273546 |
| -0.079598 -0.290320 -0.934832 -0.443210 0.758930 |
| 0.734627 -0.089070 0.079707 0.396478 0.630479 |
5 5
o22 : Matrix RR <--- RR
|
i23 : I = id_(target A)
o23 = | 1 0.000000 0.000000 0.000000 0.000000 |
| 0.000000 1 0.000000 0.000000 0.000000 |
| 0.000000 0.000000 1 0.000000 0.000000 |
| 0.000000 0.000000 0.000000 1 0.000000 |
| 0.000000 0.000000 0.000000 0.000000 1 |
5 5
o23 : Matrix RR <--- RR
|
i24 : x = solve(A,I)
o24 = | 1.374155 -0.053386 -3.245131 -0.036054 2.673772 |
| 1.220856 -1.992609 -3.776382 1.879306 1.184832 |
| -1.727056 1.673397 3.210288 -2.203965 -0.872843 |
| 0.252186 -1.421662 0.554435 1.045607 -0.750423 |
| -1.368923 0.463161 2.493175 -0.071394 -0.779715 |
5 5
o24 : Matrix RR <--- RR
|
i25 : A*x - I
o25 = 0
5 5
o25 : Matrix RR <--- RR
|
i26 : x*A - I
o26 = 0
5 5
o26 : Matrix RR <--- RR
|
Another method, which isn't generally as fast, and isn't as stable over
RR or
CC, is to lift the matrix
b along the matrix
A (see
Matrix // Matrix).
i27 : invA = I // A
o27 = | 1.374155 -0.053386 -3.245131 -0.036054 2.673772 |
| 1.220856 -1.992609 -3.776382 1.879306 1.184832 |
| -1.727056 1.673397 3.210288 -2.203965 -0.872843 |
| 0.252186 -1.421662 0.554435 1.045607 -0.750423 |
| -1.368923 0.463161 2.493175 -0.071394 -0.779715 |
5 5
o27 : Matrix RR <--- RR
|
i28 : x == invA
o28 = false
|
i29 : x - invA
o29 = | 0.000000 0.000000 0.000000 0.000000 0.000000 |
| 0.000000 0.000000 0.000000 -0.000000 0.000000 |
| 0.000000 -0.000000 -0.000000 0.000000 0.000000 |
| 0.000000 0.000000 0.000000 0.000000 0.000000 |
| 0.000000 0.000000 0.000000 0.000000 0.000000 |
5 5
o29 : Matrix RR <--- RR
|