--* From SMWATT%WATSON.vnet.ibm.com@yktvmv.watson.ibm.com  Thu Nov 10 12:34:33 1994
--* Received: from yktvmv-ob.watson.ibm.com by asharp.watson.ibm.com (AIX 3.2/UCB 5.64/930311)
--*           id AA15911; Thu, 10 Nov 1994 12:34:33 -0500
--* Received: from watson.vnet.ibm.com by yktvmv.watson.ibm.com (IBM VM SMTP V2R3)
--*    with BSMTP id 5157; Thu, 10 Nov 94 12:34:35 EST
--* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id <A.SMWATT.NOTE.VAGENT2.3649.Nov.10.12:34:21.-0500>
--*           for asbugs@watson; Thu, 10 Nov 94 12:34:33 -0500
--* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id 3623; Thu, 10 Nov 1994 12:34:21 EST
--* Received: from spadserv.watson.ibm.com by yktvmv.watson.ibm.com
--*    (IBM VM SMTP V2R3) with TCP; Thu, 10 Nov 94 12:34:15 EST
--* Received: by spadserv.watson.ibm.com (AIX 3.2/UCB 5.64/900524)
--*           id AA26707; Thu, 10 Nov 1994 12:36:39 -0500
--* Date: Thu, 10 Nov 1994 12:36:39 -0500
--* From: smwatt@spadserv.watson.ibm.com (Stephen Watt)
--* X-External-Networks: yes
--* Message-Id: <9411101736.AA26707@spadserv.watson.ibm.com>
--* To: asbugs@watson.ibm.com
--* Subject: [2] Bug: 4 constraints not checked

--@ Fixed  by:  SSD   Mon Mar 20 08:25:49 EST 1995 
--@ Tested by:  none 
--@ Summary:    Fixed type form creation for definitions with inferred types. 


-- Command line: asharp -Fo -Faso -O -v adf.as
-- Version: r1.0.2
-- Original bug file name: adf.as

--+ asharp -Fo -Faso -O -v adf.as
--+ A# (Research Version) version 1.0.2 for AIX RS/6000
--+ Bug: 4 constraints not checked
#include "aslib.as"

SI  ==> SingleInteger;

RawArray_DoubleFloat: Type == Pointer;

Array_DoubleFloat(rows: SI, cols: SI): BasicType with {
	export from DoubleFloat;

	new: (init: DoubleFloat == 0) -> %;
		++ `new(v)' creates a array with `init' in each element.

	apply: (%, SI, SI) -> DoubleFloat;
		++ `a(i,j)' extracts an element with 1-based indexing.

	set!:  (%, SI, SI, DoubleFloat) -> DoubleFloat;
		++ `a(i,j) := f' updates an element, using 1-based indexing.

	data:  % -> RawArray_DoubleFloat;
} == add {
        import {
                ArrNew: (BDFlo, BSInt) -> BArr;
                ArrElt: (BArr,  BSInt) -> BDFlo;
                ArrSet: (BArr,  BSInt, BDFlo) -> BDFlo;
        } from Builtin;

	import from DoubleFloat;

	Rep ==> BArr;

	index(i,j) ==> (rows * (i-1) + (j-1))::BSInt;

	new(init: DoubleFloat): % ==
		per ArrNew(init::BDFlo, (rows*cols)::BSInt);

	apply(a: %, i: SI, j: SI): DoubleFloat ==
		ArrElt(rep a, index(i,j))::DoubleFloat;

	set!(a: %, i: SI, j: SI, x: DoubleFloat): DoubleFloat ==
		{ ArrSet(rep a, index(i,j), x::BDFlo); x }

	-- Forget the size information associated with the type, %.
	data(a: %): RawArray_DoubleFloat ==
		a pretend RawArray_DoubleFloat;

	sample: % ==
		new(0$DoubleFloat);

	(a: %) = (b: %): Boolean == {
		for j in 1..cols repeat
			for i in 1..rows repeat
				if a(i,j) ~= b(i,j) then return false;
		true
	}

	(out: TextWriter) << (a: %): TextWriter == {
		for i in 1..rows repeat {
			out << "[";
			for j in 1..cols repeat
				out << " " << a(i,j);
			out << " ]" << newline;
		}
		out
	}
}

-- The ESSL library is in: /share/power/mathlib/essl/usr/lib

ADF ==> Array_DoubleFloat;
RAD ==> RawArray_DoubleFloat;

import from SI;

-- SVD:  A = u S vT. Updates  A := v, S := S,  B := uT B
import
	dgesvf: (
		iopt: SI,  -- "what to do" code:
			   --  0,10     singular values only
			   --  1,11     singular values and v
			   --  2,12     singular values and uT B
			   --  0, 1, 2  singular values are not ordered
		           --  10,11,12 singular values are ordered
		a:    RAD, -- input matrix
		lda:  SI,  -- leading dimension of a, as allocated
		b:    RAD, -- input matrix
		ldb:  SI,  -- leading dimension of b, as allocated
		nb:   SI,  -- number of logical columns in b
		s:    RAD, -- return vector of length n
		m:    SI,  -- conceptual row dimension of matrices a and b.
		n:    SI,  -- conceptual column dimension of matrix a
	        aux:  RAD, -- work space
		naux: SI   -- allocated dimension of aux.
			   --    naux >= 2*n + max(m, n, nb)
	) -> ()
from Foreign C "essl.h";


testSVD(): () == {
	m == 3;
	n == 3;

	a: ADF(m, n) := new 0;

	-- Invent a constructor to allow:  a:= [[...],[...],[...]]
	-- Invent another:                 a := f::ADF(3,4)
	a(1,1) := 11.0;
	a(1,2) := 12.0;
	a(1,3) := 13.0;
	a(2,1) := 21.0;
	a(2,2) := 22.0;
	a(2,3) := 23.0;
	a(3,1) := 31.0;
	a(3,2) := 32.0;
	a(3,3) := 33.0;

	nb: SI        == 1;
	b:  ADF(m,nb) := new 0;
	s:  ADF(n, 1) := new 0;

	naux == 2*n + max(n,max(m,nb));
	aux: ADF(naux, 1) := new 0;

	-- We first call the SVD routine.
	iopt := 10;
	dgesvf(iopt, data a, m, data b, m, nb, data s, m, n, data aux, naux);

	for i in 1..3 repeat
	    print << s(i,1) << newline;
}