Add sqmatrix-set routine

doc/docs/Scheme_User_Interface.md

@@ -817,7 +817,7 @@ Read/write the current eigenvectors (raw planewave amplitudes) to/from an HDF5 f
 `output-eigenvectors` is like `save-eigenvectors`, except that it saves an `evects` object returned by `get-eigenvectors`.
 Conversely `input-eigenvectors`, reads back in the eigenvectors into an `evects` object from a file that has `num-bands` bands.
 
-Currently, there's only one other interesting thing you can do with the raw eigenvectors, and that is to compute the dot-product matrix between a set of saved eigenvectors and the current eigenvectors. This can be used, e.g., to detect band crossings or to set phases consistently at different k points. The dot product is returned as a "sqmatrix" object, whose elements can be read with the `sqmatrix-size` and `sqmatrix-ref` routines.
+Currently, there's only one other interesting thing you can do with the raw eigenvectors, and that is to compute the dot-product matrix between a set of saved eigenvectors and the current eigenvectors. This can be used, e.g., to detect band crossings or to set phases consistently at different k points. The dot product is returned as a "sqmatrix" object, whose elements can be read or altered with the `sqmatrix-size`, `sqmatrix-ref`, and `sqmatrix-set` routines.
 
 **`(dot-eigenvectors ev first-band)`**  
           
@@ -831,6 +831,10 @@ Return the size *n* of an *n*x*n* sqmatrix `sm`.
           
 Return the (`i`,`j`)th element of the *n*x*n* sqmatrix *`sm`*, where {`i`,`j`} range from 0..*n*-1.
 
+**`(sqmatrix-set sm i j c)`**  
+          
+Set the (`i`,`j`)th element of the *n*x*n* sqmatrix *`sm`* to the `cnumber` *`c`*, where {`i`,`j`} range from 0..*n*-1.
+
 Inversion Symmetry
 ------------------
 

mpb/matrix-smob.c

@@ -215,6 +215,15 @@ cnumber sqmatrix_ref(SCM mo, integer i, integer j)
      return c;
 }
 
+void sqmatrix_set(SCM mo, integer i, integer j, cnumber c)
+{
+     sqmatrix *m = assert_sqmatrix_smob(mo);
+     CHECK(m && i >= 0 && j >= 0 && i < m->p && j < m->p,
+       "invalid arguments to sqmatrix-set");
+     ASSIGN_SCALAR(m->data[i * m->p + j], c.re, c.im);
+     scm_remember_upto_here_1(mo);
+}
+
 SCM sqmatrix_mult(SCM Ao, SCM Bo)
 {
     sqmatrix *A = assert_sqmatrix_smob(Ao);

mpb/mpb.scm.in

@@ -230,6 +230,8 @@
 (define-external-function sqmatrix-size false false 'integer 'SCM)
 (define-external-function sqmatrix-ref false false 'cnumber 
   'SCM 'integer 'integer)
+(define-external-function sqmatrix-set false false no-return-value
+  'SCM 'integer 'integer 'cnumber)
 (define-external-function sqmatrix-mult false false 'SCM
   'SCM 'SCM)
 (define-external-function sqmatrix-diagm false false 'SCM