c=======================================================================
c=== RBio/RBwrite mexFunction ==========================================
c=======================================================================

c RBio: a MATLAB toolbox for reading and writing sparse matrices in
c Rutherford/Boeing format.
c Copyright (c) 2006, Timothy A. Davis, Univ. Florida.  Version 1.0.


c-----------------------------------------------------------------------
c RBwrite mexFunction: write a sparse matrix to a Rutherford/Boeing file
c-----------------------------------------------------------------------
c
c   mtype = RBwrite (filename, A, Z, title, key)
c
c   A: a sparse matrix (no explicit zero entries)
c   Z: binary pattern of explicit zero entries to include in the
c	Rutherford/Boeing file.  This always has the same size as A,
c	and is always sparse.  Not used if [ ], or if nnz(Z) is 0.
c   title: title of Rutherford/Boeing file, up to 72 characters
c   key: the name of the matrix, up to 8 characters
c
c   Z is optional.  RBwrite (filename, A) uses a default 
c   title and key, and does not include any explicit zeros.
c   RBwrite (filname, A, 'title...', 'key') uses the given title and
c   key, with no Z matrix.
c
c   A must be sparse.  Z must be empty, or sparse.
c
c   mtype is a 3-character string with the file-type used
c   mtype(1):  r: 0 (real), p: 1 (pattern), c: 2 (complex),
c	       i: 3 (ineger)
c   mtype(2):  r: -1 (rectangular), u: 0 (unsymmetric), s: 1 symmetric,
c	       h: 2 (Hermitian), z: 3 (skew symmetric)
c   mtype(3):  a: assembled matrix
c-----------------------------------------------------------------------

	subroutine mexfunction (nargout, pargout, nargin, pargin)
	integer
     $	    pargout (*), pargin (*)
	integer nargout, nargin

c	----------------------------------------------------------------
c	MATLAB functions
c	----------------------------------------------------------------

	integer mxIsChar, mxClassIDFromClassName,
     $	    mxIsClass, mxIsSparse, mxIsComplex

	integer
     $	    mxGetM, mxGetN, mxGetJc, mxGetIr, mxGetPr, mxGetPi,
     $	    mxGetString, mxGetData, mxCreateNumericMatrix,
     $	    mxCreateString

c	----------------------------------------------------------------
c	local variables
c	----------------------------------------------------------------

	integer
     $	    nrow, ncol, nnz, mkind, cp, info, zmin, zmax,
     $	    skind, cmplex, wmat, cpmat, ww, kmin, kmax,
     $	    Ap, Ai, Ax, Az, Zp, Zi, Zx, Zz, znz, zrow, zcol, i,
     $	    mzkind, szkind, totcrd, ptrcrd, indcrd, valcrd,
     $	    w, valn, valn2, indn, nnz2, ptrn, i1, ititle, vals
	integer iclass
	character title*72, key*8, mtype*3, ptrfmt*20, indfmt*20,
     $	    valfmt*20, filename*1024, fmt2*20, ztype*3
	double precision t
	logical doZ, l1, is_int

c	----------------------------------------------------------------
c	check inputs
c	----------------------------------------------------------------

	if (nargin .lt. 2 .or. nargin .gt. 5. .or. nargout .gt. 2 .or.
     $	    mxIsChar (pargin (1)) .ne. 1) then
	    call mexErrMsgTxt
     $	    ('[m s] = RBwrite (filename, A, Z, title, key)')
	endif

c	----------------------------------------------------------------
c	get filename
c	----------------------------------------------------------------

	if (mxGetString (pargin (1), filename, 1024) .ne. 0) then
	    call mexErrMsgTxt ('filename too long')
	endif
	close (unit = 7)
	open (unit = 7, file = filename, err = 998)

c	----------------------------------------------------------------
c	get A
c	----------------------------------------------------------------

	if (mxIsClass (pargin (2), 'double') .ne. 1 .or.
     $	    mxIsSparse (pargin (2)) .ne. 1) then
	    call mexErrMsgTxt ('A must be sparse and double')
	endif
	cmplex = mxIsComplex (pargin (2))
	Ap = mxGetJc (pargin (2))
	Ai = mxGetIr (pargin (2))
	Ax = mxGetPr (pargin (2))
	Az = mxGetPi (pargin (2))
	nrow = mxGetM (pargin (2))
	ncol = mxGetN (pargin (2))

c	----------------------------------------------------------------
c	get title and key
c	----------------------------------------------------------------

	do 5 i = 1, 72
	    title (i:i) = ' '
5	continue
	key = '        '

	ititle = 99
	do 15 i = 3, nargin
	    if (mxIsChar (pargin (i)) .eq. 1) then
		if (ititle .eq. 99) then
c		    get the title, up to 72 characters long
		    i1 = mxGetString (pargin (i), title, 72)
		    ititle = i
		else
c		    get the key, up to 8 characters long
		    i1 = mxGetString (pargin (i), key, 8)
		endif
	    endif
15	continue

c	place a marker in the title, so we know that the
c	Rutherford/Boeing file was generated by the RBwrite mexFunction.

	title (72:72) = '|'

c	----------------------------------------------------------------
c	get Z, if present
c	----------------------------------------------------------------

	if (nargin .ge. 3 .and. ititle .gt. 3) then

	    zrow = mxGetM (pargin (3))
	    zcol = mxGetN (pargin (3))

	    if (zrow .eq. 0 .or. zcol .eq. 0) then

c		--------------------------------------------------------
c		Z matrix is empty
c		--------------------------------------------------------

		doZ = .false.

	    else

c		--------------------------------------------------------
c		get the Z matrix
c		--------------------------------------------------------

		if (mxIsClass (pargin (3), 'double') .ne. 1 .or.
     $		    mxIsSparse (pargin (3)) .ne. 1 .or.
     $		    mxIsComplex (pargin (3)) .ne. 0 .or.
     $		    zrow .ne. nrow .or. zcol .ne. ncol) then
		    call mexErrMsgTxt
     $		  ('Z must be sparse, double, real, and same size as A')
		endif
		Zp = mxGetJc (pargin (3))
		Zi = mxGetIr (pargin (3))
		Zx = mxGetPr (pargin (3))
		Zz = mxGetPi (pargin (3))
		doZ = .true.

	    endif

	else

c	    ------------------------------------------------------------
c	    no Z matrix is present
c	    ------------------------------------------------------------

	    doZ = .false.

	endif

c	----------------------------------------------------------------
c	get workspace
c	----------------------------------------------------------------

	call RBint (iclass)
	wmat  = mxCreateNumericMatrix (max (nrow, ncol)+1, 1, iclass, 0)
	cpmat = mxCreateNumericMatrix (max (nrow, ncol)+1, 1, iclass, 0)
	w = mxGetData (wmat)
	cp = mxGetData (cpmat)

c	----------------------------------------------------------------
c	determine the matrix type (RSA, RUA, etc)
c	----------------------------------------------------------------

c	find the symmetry of A (mkind, skind), and nnz(A)
	call RBkind (nrow, ncol, %val(Ap), %val(Ai), %val(Ax),
     $	    %val(Az), cmplex, mkind, skind, mtype, nnz, %val(cp),
     $	    kmin, kmax)

	if (doZ) then
c	    find the symmetry of Z and find nnz(Z)
	    call RBkind (nrow, ncol, %val(Zp), %val(Zi), %val(Zx),
     $		%val(Zz), 0, mzkind, szkind, ztype, znz, %val(cp),
     $		zmin, zmax)
	    if (znz .eq. 0) then
c		ignore the Z matrix
		doZ = .false.
	    elseif (szkind .le. 0) then
c		if Z is unsymmetric, then A+Z is unsymmetric too
		skind = szkind
	    endif
	endif

	pargout (1) = mxCreateString (mtype)

c	----------------------------------------------------------------
c	determine the required precision
c	----------------------------------------------------------------

	indfmt = '                '
	valfmt = '                    '
	is_int = mkind .eq. 3
	ww = 1
	if (mkind .ne. 1) then
	    call RBformat (nnz, %val (Ax), ww, valfmt, valn, is_int,
     $		kmin, kmax)
	    if (cmplex .eq. 1) then
		call RBformat (nnz, %val (Az), ww, valfmt, valn, is_int,
     $		kmin, kmax)
	    endif
	endif

c	----------------------------------------------------------------
c	determine the number of entries in the matrix A+Z
c	----------------------------------------------------------------

	call RBwrite (1, nrow, ncol, skind, cmplex, doZ, %val(Ap),
     $	    %val(Ai), %val(Ax), %val(Az), %val(Zp), %val(Zi), 0,
     $	    indfmt, indn, valfmt, valn, nnz2, %val(w), %val(cp))

	if (nnz2 .eq. 0) then
	    call mexErrMsgTxt ('empty matrices not handled')
	endif

c	determine pointer format.  ncol+1 integers, largest is nnz2+1
	call RBiformat (1, nnz2+1, ptrfmt, ptrn, i)
	call RBcards (ncol+1, ptrn, ptrcrd)

c	determine row index format.  nnz2 integers, largest is nrow
	call RBiformat (1, nrow, indfmt, indn, i)
	call RBcards (nnz2, indn, indcrd)

c	determine how many lines for the numerical values
	if (mkind .eq. 0 .or. mkind .eq. 3) then
c	    real or integer
	    vals = 1
	elseif (mkind .eq. 1) then
c	    pattern
	    vals = 0
	else
c	    complex
	    vals = 2
	endif
	call RBcards (vals*nnz2, valn, valcrd)

c	----------------------------------------------------------------
c	determine total number of cards
c	----------------------------------------------------------------

	totcrd = ptrcrd + indcrd + valcrd

c	----------------------------------------------------------------
c	write the header
c	----------------------------------------------------------------

        write (unit = 7, fmt = 10, err = 999)
     $          title, key,
     $          totcrd, ptrcrd, indcrd, valcrd,
     $          mtype, nrow, ncol, nnz2, 0,
     $          ptrfmt, indfmt, valfmt
10      format (a72, a8 / 4i14 / a3, 11x, 4i14 / 2a16, a20)

c	----------------------------------------------------------------
c	write the pointers
c	----------------------------------------------------------------

	call RBiflush (ptrfmt, %val (cp), ncol+1)

c	----------------------------------------------------------------
c	write the row indices
c	----------------------------------------------------------------

	call RBwrite (2, nrow, ncol, skind, cmplex, doZ, %val(Ap),
     $	    %val(Ai), %val(Ax), %val(Az), %val(Zp), %val(Zi), mkind,
     $	    indfmt, indn, valfmt, valn, nnz2, %val(w), %val(cp))

c	----------------------------------------------------------------
c	write the numerical values
c	----------------------------------------------------------------

	if (mkind .ne. 1) then
	    call RBwrite (3, nrow, ncol, skind, cmplex, doZ, %val(Ap),
     $		%val(Ai), %val(Ax), %val(Az), %val(Zp), %val(Zi), mkind,
     $		indfmt, indn, valfmt, valn, nnz2, %val(w), %val(cp))
	endif

c	----------------------------------------------------------------
c	free workspace and return
c	----------------------------------------------------------------

	close (unit = 7)
	call mxDestroyArray (%val (cpmat))
	return
998	call mexErrMsgTxt ('error openning file')
999	call mexErrMsgTxt ('error writing file')
	end