42         d=np.diagonal(ma)
43         j = np.argmax(d[k:]) + k
52         if abs(alpha) < 1.0e-18:
58         v = v.reshape((n-k-1,1))
60         ma[k+1:,k+1:] = ma[k+1:,k+1:] - np.matmul(v , np.transpose(v)) / alpha
72     d=np.diag(np.diagonal(ma))
74     return(ll,d,piv)
77 def valid(src,ll,d,piv):  argument
84     t = np.matmul(ll,np.matmul(d,np.transpose(ll)))
85     r = a - t
86     r[abs(r)<1e-10]=0.0
164 def getInvertibleMatrix(d):  argument
165   m = list(np.identity(d))
166   if d == 1:
168   if d == 2:
171     m=[[c,s],[-s,c]]
172   if d == 3:
173     m=[[0.804738, -0.310617, 0.505879], [0.505879, 
174         0.804738, -0.310617], [-0.310617, 0.505879, 0.804738]]
175   if d == 4:
178   if d == 7:
189   if d == 8:
202   if d == 9:
217   if d == 15:
254   if d == 16:
295   if d == 17:
341   if d == 32:
499   if d == 33:
669 def getDefinitePositiveMatrix(d):  argument
670     a = 1.0 * np.diag(np.array(range(1,d+1)))/d
671     p = getInvertibleMatrix(d)
674 def getSemidefinitePositiveMatrix(d,k=3):  argument
675    if d >= k + 1 :
676      a = np.diag(np.hstack([np.array(range(1,d+1-k)),np.zeros(k)])) / d
678      a = 1.0 * np.diag(np.array(range(1,d+1)))/d
679    p = getInvertibleMatrix(d)
754        r = ma - mb 
797     for d in dims:
798         ma = getInvertibleMatrix(d)
799         inp = inp + list(ma.reshape(d*d))
801         vals = vals + list(r.reshape(d*d))
815     r = np.array([0.,0.5,1.0,-0.5])
837     for d in dims:
838        ma = getDefinitePositiveMatrix(d)
839        inp = inp + list(ma.reshape(d*d))
843        vals = vals + list(l.reshape(d*d))
851        llvals = llvals + list(ll.reshape(d*d))
852        dvals = dvals + list(di.reshape(d*d))
853        permvals = permvals + list(perm.reshape(d))
855        a = np.random.randn(d*d)
857        a = a.reshape(d,d)
866        uts += list(ut.reshape(d*d)) 
867        lts += list(lt.reshape(d*d)) 
868        rndas += list(a.reshape(d*d))
870        utinvs += list(utinv.reshape(d*d)) 
871        ltinvs += list(ltinv.reshape(d*d))
872        cholinvs += list(cholinv.reshape(d*d))
904     for d in dims:
905        ma = getSemidefinitePositiveMatrix(d)
906        inp = inp + list(ma.reshape(d*d))
913        llvals = llvals + list(ll.reshape(d*d))
914        dvals = dvals + list(di.reshape(d*d))
915        permvals = permvals + list(perm.reshape(d))
1012     x = np.array([-1,0,0,0])
1033     return (np.all(np.abs(product)<1e-10))
1044     rtol = 1e-14 
1045     atol = 1e-13 
1047        rtol = 3e-2 
1048        atol = 3e-2 
1050         print(np.max(np.abs(nm-m)))
1058     eps=1e-16 
1060         eps=1e-12
1062         eps=1.0e-3
1110         #print("--------\n\n")
1118     for d in sizes:
1119         m = QR.kahan_matrix(d)
1120         thedims += [d,d,d]
1122         theMatrix += list(np.array(m).reshape(d*d))
1130         for i in range(d):
1138         theRefTau += list(np.array(tau).reshape(d))
1139         theRefR += list(np.array(r).reshape(d*d))
1140         theRefQ += list(np.array(q).reshape(d*d))
1159     p=np.identity(n)-beta * v.T .dot(v)