octave: src/DLD-FUNCTIONS/__pchip_deriv_

annotate src/DLD-FUNCTIONS/__pchip_deriv__.cc @ 7789:82be108cc558

First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code

author	David Bateman <dbateman@free.fr>
date	Sun, 27 Apr 2008 22:34:17 +0200
parents	a1dbe9d80eee
children	010e15c7be9a

rev	line source
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	1 /*
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	2
7017 a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe] jwe parents: 7016 diff changeset	3 Copyright (C) 2002, 2006, 2007 Kai Habel
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	4
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	5 This file is part of Octave.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	6
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	7 Octave is free software; you can redistribute it and/or modify it
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	8 under the terms of the GNU General Public License as published by the
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6945 diff changeset	9 Free Software Foundation; either version 3 of the License, or (at your
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6945 diff changeset	10 option) any later version.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	11
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	12 Octave is distributed in the hope that it will be useful, but WITHOUT
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	15 for more details.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	16
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	17 You should have received a copy of the GNU General Public License
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6945 diff changeset	18 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6945 diff changeset	19 <http://www.gnu.org/licenses/>.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	20
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	21 */
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	22
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	23 #ifdef HAVE_CONFIG_H
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	24 #include <config.h>
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	25 #endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	26
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	27 #include "defun-dld.h"
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	28 #include "error.h"
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	29 #include "gripes.h"
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	30 #include "oct-obj.h"
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	31 #include "utils.h"
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	32 #include "f77-fcn.h"
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	33
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	34 extern "C"
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	35 {
6242 64bad7c6a607 [project @ 2007-01-16 07:03:51 by jwe] jwe parents: 5838 diff changeset	36 F77_RET_T
64bad7c6a607 [project @ 2007-01-16 07:03:51 by jwe] jwe parents: 5838 diff changeset	37 F77_FUNC (dpchim, DPCHIM) (const octave_idx_type& n, double x, double f,
64bad7c6a607 [project @ 2007-01-16 07:03:51 by jwe] jwe parents: 5838 diff changeset	38 double *d, const octave_idx_type &incfd,
64bad7c6a607 [project @ 2007-01-16 07:03:51 by jwe] jwe parents: 5838 diff changeset	39 octave_idx_type *ierr);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	40
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	41 F77_RET_T
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	42 F77_FUNC (pchim, PCHIM) (const octave_idx_type& n, float x, float f,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	43 float *d, const octave_idx_type &incfd,
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	44 octave_idx_type *ierr);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	45 }
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	46
6945 6bbf56a9718a [project @ 2007-10-02 20:47:22 by jwe] jwe parents: 6549 diff changeset	47 // Wrapper for SLATEC/PCHIP function DPCHIM to calculate the derivates
6bbf56a9718a [project @ 2007-10-02 20:47:22 by jwe] jwe parents: 6549 diff changeset	48 // for piecewise polynomials.
6bbf56a9718a [project @ 2007-10-02 20:47:22 by jwe] jwe parents: 6549 diff changeset	49
6549 5a5a09d7deb8 [project @ 2007-04-20 06:55:29 by jwe] jwe parents: 6242 diff changeset	50 DEFUN_DLD (__pchip_deriv__, args, ,
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	51 "-- texinfo --\n\
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	52 @deftypefn {Loadable Function} {} __pchip_deriv__ (@var{x}, @var{y})\n\
6945 6bbf56a9718a [project @ 2007-10-02 20:47:22 by jwe] jwe parents: 6549 diff changeset	53 Undocumented internal function.\n\
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	54 @end deftypefn")
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	55 {
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	56 octave_value retval;
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	57 const int nargin = args.length ();
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	58
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	59 if (nargin == 2)
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	60 {
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	61 if (args(0).is_single_type () \|\| args(1).is_single_type ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	62 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	63 FloatColumnVector xvec (args(0).float_vector_value ());
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	64 FloatMatrix ymat (args(1).float_matrix_value ());
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	65
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	66 octave_idx_type nx = xvec.length ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	67 octave_idx_type nyr = ymat.rows ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	68 octave_idx_type nyc = ymat.columns ();
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	69
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	70 if (nx != nyr)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	71 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	72 error ("number of rows for x and y must match");
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	73 return retval;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	74 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	75
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	76 FloatColumnVector dvec (nx), yvec (nx);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	77 FloatMatrix dmat (nyr, nyc);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	78
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	79 octave_idx_type ierr;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	80 const octave_idx_type incfd = 1;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	81 for (int c = 0; c < nyc; c++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	82 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	83 for (int r = 0; r < nx; r++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	84 yvec(r) = ymat(r,c);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	85
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	86 F77_FUNC (pchim, PCHIM) (nx, xvec.fortran_vec (),
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	87 yvec.fortran_vec (),
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	88 dvec.fortran_vec (), incfd, &ierr);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	89
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	90 if (ierr < 0)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	91 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	92 error ("PCHIM error: %i\n", ierr);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	93 return retval;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	94 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	95
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	96 for (int r=0; r<nx; r++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	97 dmat(r,c) = dvec(r);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	98 }
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	99
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	100 retval = dmat;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	101 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	102 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	103 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	104 ColumnVector xvec (args(0).vector_value ());
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	105 Matrix ymat (args(1).matrix_value ());
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	106
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	107 octave_idx_type nx = xvec.length ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	108 octave_idx_type nyr = ymat.rows ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	109 octave_idx_type nyc = ymat.columns ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	110
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	111 if (nx != nyr)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	112 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	113 error ("number of rows for x and y must match");
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	114 return retval;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	115 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	116
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	117 ColumnVector dvec (nx), yvec (nx);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	118 Matrix dmat (nyr, nyc);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	119
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	120 octave_idx_type ierr;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	121 const octave_idx_type incfd = 1;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	122 for (int c = 0; c < nyc; c++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	123 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	124 for (int r = 0; r < nx; r++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	125 yvec(r) = ymat(r,c);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	126
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	127 F77_FUNC (dpchim, DPCHIM) (nx, xvec.fortran_vec (),
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	128 yvec.fortran_vec (),
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	129 dvec.fortran_vec (), incfd, &ierr);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	130
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	131 if (ierr < 0)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	132 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	133 error ("DPCHIM error: %i\n", ierr);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	134 return retval;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	135 }
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	136
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	137 for (int r=0; r<nx; r++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	138 dmat(r,c) = dvec(r);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	139 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	140
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	141 retval = dmat;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	142 }
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	143 }
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	144
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	145 return retval;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	146 }

Mercurial > octave

annotate src/DLD-FUNCTIONS/__pchip_deriv__.cc @ 7789:82be108cc558