From 961a8023260cbc17af40d82bd664d94b4bbb5215 Mon Sep 17 00:00:00 2001
From: Wirawan Purwanto <wirawan0@gmail.com>
Date: Tue, 24 Jan 2012 17:28:10 -0500
Subject: [PATCH] * Keep a copy of my symmetric array indexer worksheet in GIT.

  Original signature:
    -rw-r--r-- 1 wirawan0 wirawan0 1380 2010-04-27 15:55 symmetrix-array-index.py
    b873ece4610483b3cd5290c6ddbc7426  symmetrix-array-index.py
---
 math/symmetrix-array-index.PY | 73 +++++++++++++++++++++++++++++++++++
 1 file changed, 73 insertions(+)
 create mode 100644 math/symmetrix-array-index.PY

diff --git a/math/symmetrix-array-index.PY b/math/symmetrix-array-index.PY
new file mode 100644
index 0000000..54da7ba
--- /dev/null
+++ b/math/symmetrix-array-index.PY
@@ -0,0 +1,73 @@
+# 20100427
+
+"""
+Symmetric array layout (lower diagonal, Fortran column-major ordering):
+
+
+     ----> j
+  |  (1,1)
+i |  (2,1)   (2,2)
+  v  (3,1)   (3,2)   (3,3)
+       :       :       :
+     (N,1)   (N,2)   (N,3) ... (N,N)
+
+In linear form:
+
+       1
+       2      N+1
+       3      N+2     2N
+       :       :       :
+       N     N+N-1
+
+
+
+iskip = i - j + 1   (easy)
+
+jskip determination
+  i     jskip        jskip
+        (cumulant)   (total)
+  1        0
+  2        N
+  3        N-1
+  4        N-2
+  ...
+  j        N-(j-2)   (j-1)N - (j-2)(j-1)/2
+  ...
+
+
+Note:
+
+    {m}
+    SUM a = 1 + 2 + ... + m = m(m+1)/2
+   {a=1}
+"""
+
+import numpy
+
+def test_jskip_1(N, M):
+  jskip1 = 0
+  jskip2 = 0
+  for m in xrange(1, M+1):
+    if m == 1:
+      cum = 0
+    else:
+      cum = N - (m-2)
+    jskip2 = (m-1)*N - (m-2)*(m-1)/2
+    jskip1 += cum
+    print "%5d %5d %5d %5d" % (m, cum, jskip1, jskip2)
+
+
+def copy_over_array(N, arr_L):
+  rslt = numpy.zeros((N,N))
+  for i in xrange(N):
+    for j in xrange(N):
+      if j > i: continue
+      ii = i+1 # Fortranize it
+      jj = j+1 # Fortranize it
+      jskip = (jj-1)*N - (jj-2)*(jj-1)/2
+      iskip = ii - jj + 1
+      ldiag_index = iskip + jskip - 1 # -1 for C-izing again
+      # indices printed in Fortran 1-based indices
+      print "%5d %5d %5d %5d %5d" % (ii,jj,iskip,jskip,ldiag_index+1)
+      rslt[i,j] = arr_L[ldiag_index]
+  return rslt