Add (real&complex) Schur decomposition

configure

@@ -5389,7 +5389,7 @@ fi
 
 
     vars="vectclapi.c vectcl.c linalg.c arrayshape.c nacomplex.c hsfft.c fft.c
-				svd.c eig.c
+				svd.c eig.c schur.c
 				bcexecute.c vmparser.c
 				clapack_cutdown.c tcl_xerbla.c"
     for i in $vars; do

configure.in

@@ -75,7 +75,7 @@ AC_TYPE_INTPTR_T
 #-----------------------------------------------------------------------
 
 TEA_ADD_SOURCES([vectclapi.c vectcl.c linalg.c arrayshape.c nacomplex.c hsfft.c fft.c 
-				svd.c eig.c
+				svd.c eig.c schur.c
 				bcexecute.c vmparser.c 
 				clapack_cutdown.c tcl_xerbla.c])
 TEA_ADD_HEADERS([generic/vectcl.h generic/nacomplex.h generic/hsfft.h])

generic/clapack_cutdown.h

@@ -16,4 +16,6 @@ MODULE_SCOPE /* Subroutine */ int dgesv_ (Tcl_Interp *interp, integer *n, intege
 MODULE_SCOPE /* Subroutine */ int zgesv_ (Tcl_Interp *interp, integer *n, integer *nrhs, doublecomplex *a, 	integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer *	info);
 MODULE_SCOPE /* Subroutine */ int dgesvx_ (Tcl_Interp *interp, char *fact, char *trans, integer *n, integer *	nrhs, doublereal *a, integer *lda, doublereal *af, integer *ldaf, 	integer *ipiv, char *equed, doublereal *r__, doublereal *c__, 	doublereal *b, integer *ldb, doublereal *x, integer *ldx, doublereal *	rcond, doublereal *ferr, doublereal *berr, doublereal *work, integer *	iwork, integer *info);
 MODULE_SCOPE /* Subroutine */ int zgesvx_ (Tcl_Interp *interp, char *fact, char *trans, integer *n, integer *	nrhs, doublecomplex *a, integer *lda, doublecomplex *af, integer *	ldaf, integer *ipiv, char *equed, doublereal *r__, doublereal *c__, 	doublecomplex *b, integer *ldb, doublecomplex *x, integer *ldx, 	doublereal *rcond, doublereal *ferr, doublereal *berr, doublecomplex *	work, doublereal *rwork, integer *info);
+MODULE_SCOPE /* Subroutine */ int dgees_ (Tcl_Interp *interp, char *jobvs, char *sort, L_fp select, integer *n, 	doublereal *a, integer *lda, integer *sdim, doublereal *wr, 	doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work, 	integer *lwork, logical *bwork, integer *info);
+MODULE_SCOPE /* Subroutine */ int zgees_ (Tcl_Interp *interp, char *jobvs, char *sort, L_fp select, integer *n, 	doublecomplex *a, integer *lda, integer *sdim, doublecomplex *w, 	doublecomplex *vs, integer *ldvs, doublecomplex *work, integer *lwork, 	 doublereal *rwork, logical *bwork, integer *info);
 

generic/nacomplex.c

@@ -2,8 +2,6 @@
 #include <stdlib.h>
 #include <string.h>
 
-/* TODO: These are just stubs to test the different types */
-
 /*
  * Functions hndling the Tcl_ObjType Complex
  */

generic/schur.c

@@ -0,0 +1,182 @@
+#include "schur.h"
+#include "clapack_cutdown.h"
+
+#define MIN(X, Y) ((X)<(Y) ? X : Y)
+#define MAX(X, Y) ((X)>(Y) ? X : Y)
+
+static int doSchur(Tcl_Interp *interp, Tcl_Obj *matrix, Tcl_Obj **Z,  Tcl_Obj **T) {
+	/* Compute Schur decomposition of a matrix.
+	 * Return Schur vectors in Z and Schur form in T, 
+	 */
+
+	/* Convert matrix to VecTcl object */
+	NumArrayInfo *info = NumArrayGetInfoFromObj(interp, matrix);
+	if (!info) { return TCL_ERROR; }
+
+	/* Check that it is a square matrix */
+	if (info->nDim != 2) {
+		/* Could be a scalar. In this case return the trivial
+		 * decomposition */
+		if (ISSCALARINFO(info)) {
+			*T = Tcl_DuplicateObj(matrix);
+			*Z = Tcl_NewDoubleObj(1.0);
+			return TCL_OK;
+		}
+
+		Tcl_SetResult(interp, "Schur decomposition is only defined for square matrix", NULL);
+		return TCL_ERROR;
+	}
+
+
+	/* get matrix dimensions */
+	long int m = info->dims[0];
+	long int n = info->dims[1];
+
+	if (m != n) {
+		Tcl_SetResult(interp, "Schur decomposition is only defined for square matrix", NULL);
+		return TCL_ERROR;
+	}
+
+	char *jobvs = "V";
+	char *sort = "N";
+
+	if (info->type != NumArray_Complex128) {
+		/* Real-valued matrix, prepare for dgees */
+		/* create a column-major copy of matrix 
+		 * This also converts an integer matrix to double */
+		*T = NumArrayNewMatrixColMaj(NumArray_Float64, m, n);
+		NumArrayObjCopy(interp, matrix, *T);
+
+		*Z = NumArrayNewMatrixColMaj(NumArray_Float64, m, m);
+
+		/* Extract the raw pointers from the VecTcl objects */
+		double *Tptr = NumArrayGetPtrFromObj(interp, *T);
+		double *Zptr = NumArrayGetPtrFromObj(interp, *Z);
+
+		/* Space to store the eigenvalues */
+		doublereal *wr = ckalloc(sizeof(doublereal)*n);
+		doublereal *wi = ckalloc(sizeof(doublereal)*n);
+
+		/* setup workspace arrays */
+		integer lwork = 3*n;
+		doublereal* work=ckalloc(sizeof(doublereal)*lwork);
+		logical *bwork = NULL;
+		integer sdim=0;
+
+		/* Leading dimensions of T and Vr 
+		 * Don't compute left vectors. */
+		integer ldt = n;
+		integer ldz = n;
+		integer info;
+
+/* Subroutine  int dgees_(char *jobvs, char *sort, L_fp select, 
+ * integer *n, doublereal *a, integer *lda, integer *sdim, 
+ * doublereal *wr, doublereal *wi, doublereal *vs, integer *ldvs, 
+ * doublereal *work, integer *lwork, logical *bwork, integer *info) */
+
+
+		/* call out to dgees */
+		int errcode=dgees_(interp, jobvs, sort, NULL, 
+				&n, Tptr, &ldt, &sdim,
+				wr, wi, Zptr, &ldz,
+				work, &lwork, bwork, &info);
+
+		/* free workspace */
+		ckfree(work);
+		ckfree(wr); ckfree(wi);
+
+		if (errcode != TCL_OK) {
+			/* release temporary storage for result */
+			Tcl_DecrRefCount(*Z);
+			Tcl_DecrRefCount(*T);
+			if (errcode > 0) {
+				RESULTPRINTF(("DGEES failed to converge at eigenvector %d ", info));
+			}
+			return TCL_ERROR;
+		}
+
+		return TCL_OK;
+
+	} else {
+		/* Complex matrix, prepare for zgees */
+		/* create a column-major copy of matrix 
+		 * This also converts an integer matrix to double */
+		*T = NumArrayNewMatrixColMaj(NumArray_Complex128, m, n);
+		NumArrayObjCopy(interp, matrix, *T);
+
+		*Z = NumArrayNewMatrixColMaj(NumArray_Complex128, m, m);
+
+		/* Extract the raw pointers from the VecTcl objects */
+		doublecomplex *Tptr = NumArrayGetPtrFromObj(interp, *T);
+		doublecomplex *Zptr = NumArrayGetPtrFromObj(interp, *Z);
+
+		/* Space to store the eigenvalues */
+		doublecomplex *w = ckalloc(sizeof(doublecomplex)*n);
+
+		/* setup workspace arrays */
+		integer lwork = 2*n;
+		doublecomplex *work=ckalloc(sizeof(doublecomplex)*lwork);
+		doublereal *rwork=ckalloc(sizeof(doublereal)*n);
+		logical *bwork = NULL;
+		integer sdim=0;
+
+		/* Leading dimensions of T and Vr 
+		 * Don't compute left vectors. */
+		integer ldt = n;
+		integer ldz = n;
+		integer info;
+
+	/* Subroutine int zgees_(char *jobvs, char *sort, L_fp select, 
+	 * integer *n, doublecomplex *a, integer *lda, integer *sdim, 
+	 * doublecomplex *w, doublecomplex *vs, integer *ldvs, 
+	 * doublecomplex *work, integer *lwork, doublereal *rwork, logical *bwork, integer *info) */
+
+		/* call out to dgees */
+		int errcode=zgees_(interp, jobvs, sort, NULL, 
+				&n, Tptr, &ldt, &sdim,
+				w, Zptr, &ldz,
+				work, &lwork, rwork, bwork, &info);
+
+		/* free workspace */
+		ckfree(work);
+		ckfree(rwork);
+		ckfree(w);
+
+		if (errcode != TCL_OK) {
+			/* release temporary storage for result */
+			Tcl_DecrRefCount(*Z);
+			Tcl_DecrRefCount(*T);
+			if (errcode > 0) {
+				RESULTPRINTF(("ZGEES failed to converge at eigenvector %d ", info));
+			}
+			return TCL_ERROR;
+		}
+
+		return TCL_OK;
+
+	}
+}
+
+int NumArraySchurCmd(ClientData dummy, Tcl_Interp *interp, int objc, Tcl_Obj *const *objv) {
+	if (objc != 2) {
+		Tcl_WrongNumArgs(interp, 1, objv, "matrix");
+		return TCL_ERROR;
+	}
+
+	Tcl_Obj *matrix = objv[1];
+
+	Tcl_Obj *Z, *T;
+
+	if (doSchur(interp, matrix, &Z, &T) != TCL_OK) {
+		return TCL_ERROR;
+	}
+
+	/* return as list */
+	Tcl_Obj *result=Tcl_NewObj();
+	Tcl_ListObjAppendElement(interp, result, Z);
+	Tcl_ListObjAppendElement(interp, result, T);
+
+	Tcl_SetObjResult(interp, result);
+	return TCL_OK;
+}
+

generic/schur.h

@@ -0,0 +1,12 @@
+/* function definitions for basic linear algebra
+ * matrix decompositions / equation system solving */
+
+#ifndef SCHUR_H
+#define SCHUR_H
+#include "vectclInt.h"
+
+/* Compute the Schur form */
+int NumArraySchurCmd(ClientData dummy, Tcl_Interp *interp,
+    int objc, Tcl_Obj *const *objv);
+
+#endif

generic/vectcl.c

@@ -8,6 +8,7 @@
 #include "fft.h"
 #include "svd.h"
 #include "eig.h"
+#include "schur.h"
 #include "bcexecute.h"
 #include "vmparser.h"
 
@@ -245,6 +246,7 @@ static const EnsembleMap implementationMap[] = {
 	{"eig", NumArrayEigCmd, NULL},
 	{"svd1", NumArraySVD1Cmd, NULL},
 	{"svd", NumArraySVDCmd, NULL},
+	{"schur", NumArraySchurCmd, NULL},
 	/* Reductions */
 	{"sum", NumArraySumCmd, NULL},
 	{"axismin", NumArrayAxisMinCmd, NULL},

tools/cherrypick_clapack.tcl

@@ -4,8 +4,8 @@
 
 set CLAPACK_DIR /Users/chris/Sources/CLAPACK-3.2.1/
 
-set tooldir [file dirname [info script]]
-set destfile [file join $tooldir ../../generic/clackap_cutdown.c]
+set tooldir [file normalize [file dirname [info script]]]
+set destfile [file normalize [file join $tooldir ../generic/clapack_cutdown.c]]
 
 # helper procs for parsing
 proc shift {} {
@@ -422,5 +422,6 @@ proc run_generator {} {
 	dgesdd_ zgesdd_ dgemm_ zgemm_ \
 	dsyevr_ zheevr_ dgeev_ zgeev_ \
 	dgelss_ zgelss_ dgelsy_ zgelsy_ \
-	dgesv_ zgesv_ dgesvx_ zgesvx_
+	dgesv_ zgesv_ dgesvx_ zgesvx_ \
+	dgees_ zgees_ 
 }