Merged in implement_more_math_mod (pull request #18)

Implement more of the math module

Aaron R Robinson committed 9 years ago

Commit a000a972 parents be65a2f 9e613375 files · +792 −15Browse files

Unified · Split

5 files

readme.md +51 −0
src/vm/disvm.gyp +2 −0
src/vm/math/Mathmod.cpp +215 −15
src/vm/math/dgemm.c +411 −0
src/vm/math/lsame.c +113 −0

@@ -116,6 +116,57 @@ The above copyright notice and this permission notice shall be included in all c · 115 unmodified lines
116	116
117	117	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
118	118
	119	+### [LAPACK](http://www.netlib.org/lapack/#_licensing) ###
	120	+
	121	+Copyright (c) 1992-2013 The University of Tennessee and The University
	122	+ of Tennessee Research Foundation. All rights
	123	+ reserved.
	124	+Copyright (c) 2000-2013 The University of California Berkeley. All
	125	+ rights reserved.
	126	+Copyright (c) 2006-2013 The University of Colorado Denver. All rights
	127	+ reserved.
	128	+
	129	+$COPYRIGHT$
	130	+
	131	+Additional copyrights may follow
	132	+
	133	+$HEADER$
	134	+
	135	+Redistribution and use in source and binary forms, with or without
	136	+modification, are permitted provided that the following conditions are
	137	+met:
	138	+
	139	+- Redistributions of source code must retain the above copyright
	140	+ notice, this list of conditions and the following disclaimer.
	141	+
	142	+- Redistributions in binary form must reproduce the above copyright
	143	+ notice, this list of conditions and the following disclaimer listed
	144	+ in this license in the documentation and/or other materials
	145	+ provided with the distribution.
	146	+
	147	+- Neither the name of the copyright holders nor the names of its
	148	+ contributors may be used to endorse or promote products derived from
	149	+ this software without specific prior written permission.
	150	+
	151	+The copyright holders provide no reassurances that the source code
	152	+provided does not infringe any patent, copyright, or any other
	153	+intellectual property rights of third parties. The copyright holders
	154	+disclaim any liability to any recipient for claims brought against
	155	+recipient by any third party for infringement of that parties
	156	+intellectual property rights.
	157	+
	158	+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
	159	+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
	160	+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
	161	+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
	162	+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
	163	+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
	164	+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
	165	+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
	166	+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
	167	+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
	168	+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
	169	+
119	170	### Unless otherwise stated in a `license.txt` or above - [MIT License](https://opensource.org/licenses/MIT) ###
120	171
121	172	Copyright (c) 2017 arr

src/vm/disvm.gyp+2−0

@@ -24,6 +24,8 @@ · 23 unmodified lines
24	24	'vm_memory.cpp',
25	25	'vm_exception_handler.cpp',
26	26	'math/Mathmod.cpp',
	27	+ 'math/dgemm.c',
	28	+ 'math/lsame.c',
27	29	'sys/print.cpp',
28	30	'sys/Sysmod.cpp',
29	31	'sys/fd_manager.cpp',

src/vm/math/Mathmod.cpp+215−15

@@ -14,6 +14,7 @@ using disvm::debug::log_level_t; · 13 unmodified lines
14	14
15	15	using disvm::runtime::dereference_nil;
16	16	using disvm::runtime::intrinsic_type_desc;
	17	+using disvm::runtime::out_of_range_memory;
17	18	using disvm::runtime::real_t;
18	19	using disvm::runtime::vm_alloc_t;
19	20	using disvm::runtime::vm_array_t;
@@ -119,6 +120,8 @@ void · 99 unmodified lines
119	120	Math_FPstatus(vm_registers_t &r, vm_t &vm)
120	121	{
121	122	auto &fp = r.stack.peek_frame()->base<F_Math_FPstatus>();
	123	+
	124	+ // [PAL] Properly implement floating point control
122	125	throw vm_system_exception{ "Function not implemented" };
123	126	}
124	127
@@ -290,7 +293,28 @@ void · 165 unmodified lines
290	293	Math_export_int(vm_registers_t &r, vm_t &vm)
291	294	{
292	295	auto &fp = r.stack.peek_frame()->base<F_Math_export_int>();
293		− throw vm_system_exception{ "Function not implemented" };
	296	+
	297	+ auto buffer = vm_alloc_t::from_allocation<vm_array_t>(fp.b);
	298	+ assert(buffer->get_element_type()->is_equal(intrinsic_type_desc::type<byte_t>().get()));
	299	+
	300	+ auto values = vm_alloc_t::from_allocation<vm_array_t>(fp.x);
	301	+ assert(values->get_element_type()->is_equal(intrinsic_type_desc::type<word_t>().get()));
	302	+ const auto values_len = values->get_length();
	303	+
	304	+ if (buffer->get_length() != (sizeof(word_t) * values_len))
	305	+ throw vm_user_exception{ "Invalid buffer size based on conversion" };
	306	+
	307	+ // Insert as big endian
	308	+ for (int i = 0; i < values_len; ++i)
	309	+ {
	310	+ auto x = values->at<word_t>(i);
	311	+
	312	+ auto base = reinterpret_cast<byte_t >(buffer->at(i sizeof(x)));
	313	+ base[0] = static_cast<byte_t>(x >> 24);
	314	+ base[1] = static_cast<byte_t>(x >> 16);
	315	+ base[2] = static_cast<byte_t>(x >> 8);
	316	+ base[3] = static_cast<byte_t>(x);
	317	+ }
294	318	}
295	319
296	320	void
@@ -305,7 +329,7 @@ Math_export_real(vm_registers_t &r, vm_t &vm) · 8 unmodified lines
305	329	assert(values->get_element_type()->is_equal(intrinsic_type_desc::type<real_t>().get()));
306	330	const auto values_len = values->get_length();
307	331
308		− if (buffer->get_length() != (8 * values_len))
	332	+ if (buffer->get_length() != (sizeof(real_t) * values_len))
309	333	throw vm_user_exception{ "Invalid buffer size based on conversion" };
310	334
311	335	// Insert as big endian
@@ -332,7 +356,32 @@ void · 20 unmodified lines
332	356	Math_export_real32(vm_registers_t &r, vm_t &vm)
333	357	{
334	358	auto &fp = r.stack.peek_frame()->base<F_Math_export_real32>();
335		− throw vm_system_exception{ "Function not implemented" };
	359	+
	360	+ auto buffer = vm_alloc_t::from_allocation<vm_array_t>(fp.b);
	361	+ assert(buffer->get_element_type()->is_equal(intrinsic_type_desc::type<byte_t>().get()));
	362	+
	363	+ auto values = vm_alloc_t::from_allocation<vm_array_t>(fp.x);
	364	+ assert(values->get_element_type()->is_equal(intrinsic_type_desc::type<real_t>().get()));
	365	+ const auto values_len = values->get_length();
	366	+
	367	+ if (buffer->get_length() != (sizeof(float) * values_len))
	368	+ throw vm_user_exception{ "Invalid buffer size based on conversion" };
	369	+
	370	+ // Insert as big endian
	371	+ for (int i = 0; i < values_len; ++i)
	372	+ {
	373	+ auto x_r = values->at<real_t>(i);
	374	+ auto x = static_cast<float>(x_r);
	375	+ uint32_t t;
	376	+ std::memcpy(&t, &x, sizeof(x));
	377	+ static_assert(sizeof(t) == sizeof(x), "Short-real should be 32-bits");
	378	+
	379	+ auto base = reinterpret_cast<byte_t >(buffer->at(i sizeof(x)));
	380	+ base[0] = static_cast<byte_t>(t >> 24);
	381	+ base[1] = static_cast<byte_t>(t >> 16);
	382	+ base[2] = static_cast<byte_t>(t >> 8);
	383	+ base[3] = static_cast<byte_t>(t);
	384	+ }
336	385	}
337	386
338	387	void
@@ -384,17 +433,79 @@ Math_fmod(vm_registers_t &r, vm_t &vm) · 45 unmodified lines
384	433	*fp.ret = std::fmod(fp.x, fp.y);
385	434	}
386	435
	436	+// Forward declaration
	437	+extern "C" int dgemm_(
	438	+ char transa, char transb,
	439	+ int m, int n, int *k,
	440	+ double *alpha,
	441	+ double a, int lda,
	442	+ double b, int ldb,
	443	+ double *beta,
	444	+ double c, int ldc,
	445	+ int a_len, int b_len, int c_len);
	446	+
387	447	void
388	448	Math_gemm(vm_registers_t &r, vm_t &vm)
389	449	{
390	450	auto &fp = r.stack.peek_frame()->base<F_Math_gemm>();
391		− throw vm_system_exception{ "Function not implemented" };
	451	+
	452	+ char t_a[2] = {};
	453	+ t_a[0] = static_cast<char>(fp.transa);
	454	+
	455	+ char t_b[2] = {};
	456	+ t_b[0] = static_cast<char>(fp.transb);
	457	+
	458	+ auto m = fp.m;
	459	+ auto n = fp.n;
	460	+ auto k = fp.k;
	461	+
	462	+ auto alpha = fp.alpha;
	463	+ auto beta = fp.beta;
	464	+
	465	+ auto lda = fp.lda;
	466	+ auto a_v = vm_alloc_t::from_allocation<vm_array_t>(fp.a);
	467	+ if (a_v == nullptr)
	468	+ throw dereference_nil{ "gemm - a" };
	469	+
	470	+ auto ldb = fp.ldb;
	471	+ auto b_v = vm_alloc_t::from_allocation<vm_array_t>(fp.b);
	472	+ if (b_v == nullptr)
	473	+ throw dereference_nil{ "gemm - b" };
	474	+
	475	+ auto ldc = fp.ldc;
	476	+ auto c_v = vm_alloc_t::from_allocation<vm_array_t>(fp.c);
	477	+ if (c_v == nullptr)
	478	+ throw dereference_nil{ "gemm - c" };
	479	+
	480	+ assert(a_v->get_element_type()->is_equal(intrinsic_type_desc::type<real_t>().get()));
	481	+ assert(b_v->get_element_type()->is_equal(intrinsic_type_desc::type<real_t>().get()));
	482	+ assert(c_v->get_element_type()->is_equal(intrinsic_type_desc::type<real_t>().get()));
	483	+
	484	+ const auto a_len = a_v->get_length();
	485	+ const auto b_len = b_v->get_length();
	486	+ const auto c_len = c_v->get_length();
	487	+
	488	+ // Call BLAS version of function
	489	+ const auto result = dgemm_(
	490	+ t_a, t_b,
	491	+ &m, &n, &k,
	492	+ &alpha,
	493	+ reinterpret_cast<double *>(a_v->at(0)), &lda,
	494	+ reinterpret_cast<double *>(b_v->at(0)), &ldb,
	495	+ &beta,
	496	+ reinterpret_cast<double *>(c_v->at(0)), &ldc,
	497	+ a_len, b_len, c_len);
	498	+
	499	+ if (result != 0)
	500	+ throw vm_user_exception{ "gemm invalid argument" };
392	501	}
393	502
394	503	void
395	504	Math_getFPcontrol(vm_registers_t &r, vm_t &vm)
396	505	{
397	506	auto &fp = r.stack.peek_frame()->base<F_Math_getFPcontrol>();
	507	+
	508	+ // [PAL] Properly implement floating point control
398	509	throw vm_system_exception{ "Function not implemented" };
399	510	}
400	511
@@ -402,6 +513,8 @@ void · 1 unmodified line
402	513	Math_getFPstatus(vm_registers_t &r, vm_t &vm)
403	514	{
404	515	auto &fp = r.stack.peek_frame()->base<F_Math_getFPstatus>();
	516	+
	517	+ // [PAL] Properly implement floating point control
405	518	throw vm_system_exception{ "Function not implemented" };
406	519	}
407	520
@@ -416,7 +529,34 @@ void · 8 unmodified lines
416	529	Math_iamax(vm_registers_t &r, vm_t &vm)
417	530	{
418	531	auto &fp = r.stack.peek_frame()->base<F_Math_iamax>();
419		− throw vm_system_exception{ "Function not implemented" };
	532	+
	533	+ // [SPEC] According to the Inferno implementation (libmath/blas.c), 0 is returned
	534	+ // if the array is empty. This means there is no indication the supplied array was
	535	+ // empty or the first index was max. The spec does not describe this condition,
	536	+ // therefore this implementation will indicate with -1.
	537	+ *fp.ret = -1;
	538	+
	539	+ auto result = vm_alloc_t::from_allocation<vm_array_t>(fp.x);
	540	+ assert(result->get_element_type()->is_equal(intrinsic_type_desc::type<real_t>().get()));
	541	+ if (result->get_length() <= 0)
	542	+ return;
	543	+
	544	+ auto r_raw = reinterpret_cast<real_t *>(result->at(0));
	545	+
	546	+ auto max_idx = 0;
	547	+ auto r_len = result->get_length();
	548	+ auto curr_max = std::fabs(r_raw[max_idx]);
	549	+ for (auto i = word_t{ 1 }; i < r_len; ++i)
	550	+ {
	551	+ auto t = std::fabs(r_raw[i]);
	552	+ if (curr_max < t)
	553	+ {
	554	+ max_idx = i;
	555	+ curr_max = t;
	556	+ }
	557	+ }
	558	+
	559	+ *fp.ret = max_idx;
420	560	}
421	561
422	562	void
@@ -430,7 +570,27 @@ void · 7 unmodified lines
430	570	Math_import_int(vm_registers_t &r, vm_t &vm)
431	571	{
432	572	auto &fp = r.stack.peek_frame()->base<F_Math_import_int>();
433		− throw vm_system_exception{ "Function not implemented" };
	573	+
	574	+ auto buffer = vm_alloc_t::from_allocation<vm_array_t>(fp.b);
	575	+ assert(buffer->get_element_type()->is_equal(intrinsic_type_desc::type<byte_t>().get()));
	576	+
	577	+ auto result = vm_alloc_t::from_allocation<vm_array_t>(fp.x);
	578	+ assert(result->get_element_type()->is_equal(intrinsic_type_desc::type<word_t>().get()));
	579	+ const auto result_len = result->get_length();
	580	+
	581	+ if (buffer->get_length() != (sizeof(word_t) * result_len))
	582	+ throw vm_user_exception{ "Invalid buffer size based on conversion" };
	583	+
	584	+ // Process as big endian
	585	+ for (int i = 0; i < result_len; ++i)
	586	+ {
	587	+ uint32_t t = buffer->at<byte_t>(0);
	588	+ for (int j = 1; j < sizeof(word_t); ++j)
	589	+ t = (t << 8) \| buffer->at<byte_t>(j);
	590	+
	591	+ auto x = static_cast<word_t>(t);
	592	+ result->at<word_t>(i) = x;
	593	+ }
434	594	}
435	595
436	596	void
@@ -445,14 +605,14 @@ Math_import_real(vm_registers_t &r, vm_t &vm) · 8 unmodified lines
445	605	assert(result->get_element_type()->is_equal(intrinsic_type_desc::type<real_t>().get()));
446	606	const auto result_len = result->get_length();
447	607
448		− if (buffer->get_length() != (8 * result_len))
	608	+ if (buffer->get_length() != (sizeof(real_t) * result_len))
449	609	throw vm_user_exception{ "Invalid buffer size based on conversion" };
450	610
451	611	// Process as big endian
452	612	for (int i = 0; i < result_len; ++i)
453	613	{
454	614	uint64_t t = buffer->at<byte_t>(0);
455		− for (int j = 1; j < 8; ++j)
	615	+ for (int j = 1; j < sizeof(real_t); ++j)
456	616	t = (t << 8) \| buffer->at<byte_t>(j);
457	617
458	618	real_t x;
@@ -467,7 +627,31 @@ void · 8 unmodified lines
467	627	Math_import_real32(vm_registers_t &r, vm_t &vm)
468	628	{
469	629	auto &fp = r.stack.peek_frame()->base<F_Math_import_real32>();
470		− throw vm_system_exception{ "Function not implemented" };
	630	+
	631	+ auto buffer = vm_alloc_t::from_allocation<vm_array_t>(fp.b);
	632	+ assert(buffer->get_element_type()->is_equal(intrinsic_type_desc::type<byte_t>().get()));
	633	+
	634	+ auto result = vm_alloc_t::from_allocation<vm_array_t>(fp.x);
	635	+ assert(result->get_element_type()->is_equal(intrinsic_type_desc::type<real_t>().get()));
	636	+ const auto result_len = result->get_length();
	637	+
	638	+ if (buffer->get_length() != (sizeof(float) * result_len))
	639	+ throw vm_user_exception{ "Invalid buffer size based on conversion" };
	640	+
	641	+ // Process as big endian
	642	+ for (int i = 0; i < result_len; ++i)
	643	+ {
	644	+ uint32_t t = buffer->at<byte_t>(0);
	645	+ for (int j = 1; j < sizeof(float); ++j)
	646	+ t = (t << 8) \| buffer->at<byte_t>(j);
	647	+
	648	+ float x_f;
	649	+ std::memcpy(&x_f, &t, sizeof(t));
	650	+ static_assert(sizeof(t) == sizeof(x_f), "Short-real should be 32-bits");
	651	+
	652	+ auto x = real_t{ x_f };
	653	+ result->at<real_t>(i) = x;
	654	+ }
471	655	}
472	656
473	657	void
@@ -481,6 +665,8 @@ void · 7 unmodified lines
481	665	Math_j0(vm_registers_t &r, vm_t &vm)
482	666	{
483	667	auto &fp = r.stack.peek_frame()->base<F_Math_j0>();
	668	+
	669	+ // [PAL] Consume platform API
484	670	throw vm_system_exception{ "Function not implemented" };
485	671	}
486	672
@@ -488,6 +674,8 @@ void · 1 unmodified line
488	674	Math_j1(vm_registers_t &r, vm_t &vm)
489	675	{
490	676	auto &fp = r.stack.peek_frame()->base<F_Math_j1>();
	677	+
	678	+ // [PAL] Consume platform API
491	679	throw vm_system_exception{ "Function not implemented" };
492	680	}
493	681
@@ -495,6 +683,8 @@ void · 1 unmodified line
495	683	Math_jn(vm_registers_t &r, vm_t &vm)
496	684	{
497	685	auto &fp = r.stack.peek_frame()->base<F_Math_jn>();
	686	+
	687	+ // [PAL] Consume platform API
498	688	throw vm_system_exception{ "Function not implemented" };
499	689	}
500	690
@@ -530,7 +720,13 @@ void · 29 unmodified lines
530	720	Math_modf(vm_registers_t &r, vm_t &vm)
531	721	{
532	722	auto &fp = r.stack.peek_frame()->base<F_Math_modf>();
533		− throw vm_system_exception{ "Function not implemented" };
	723	+
	724	+ real_t integral_out;
	725	+ const auto fractional_part = std::modf(fp.x, &integral_out);
	726	+ const auto integral_part = static_cast<word_t>(integral_out);
	727	+
	728	+ fp.ret->t0 = integral_part;
	729	+ fp.ret->t1 = fractional_part;
534	730	}
535	731
536	732	void
@@ -565,9 +761,7 @@ Math_norm1(vm_registers_t &r, vm_t &vm) · 28 unmodified lines
565	761	void
566	762	Math_norm2(vm_registers_t &r, vm_t &vm)
567	763	{
568		− // [SPEC] This function is defined as: norm2(x) = sqrt(dot(x, x))
569		− // The Inferno implementation (libmath/blas.c) doesn't appear to take
570		− // the square root of the resulting dot product.
	764	+ // Defined as: norm2(x) = sqrt(dot(x, x))
571	765	auto &fp = r.stack.peek_frame()->base<F_Math_norm2>();
572	766
573	767	auto xs = vm_alloc_t::from_allocation<vm_array_t>(fp.x);
@@ -600,7 +794,7 @@ void · 26 unmodified lines
600	794	Math_pow10(vm_registers_t &r, vm_t &vm)
601	795	{
602	796	auto &fp = r.stack.peek_frame()->base<F_Math_pow10>();
603		− throw vm_system_exception{ "Function not implemented" };
	797	+ *fp.ret = std::pow(10., static_cast<real_t>(fp.p));
604	798	}
605	799
606	800	void
@@ -646,7 +840,7 @@ void · 39 unmodified lines
646	840	Math_scalbn(vm_registers_t &r, vm_t &vm)
647	841	{
648	842	auto &fp = r.stack.peek_frame()->base<F_Math_scalbn>();
649		− throw vm_system_exception{ "Function not implemented" };
	843	+ *fp.ret = std::scalbln(fp.x, fp.n);
650	844	}
651	845
652	846	void
@@ -723,6 +917,8 @@ void · 70 unmodified lines
723	917	Math_y0(vm_registers_t &r, vm_t &vm)
724	918	{
725	919	auto &fp = r.stack.peek_frame()->base<F_Math_y0>();
	920	+
	921	+ // [PAL] Consume platform API
726	922	throw vm_system_exception{ "Function not implemented" };
727	923	}
728	924
@@ -730,6 +926,8 @@ void · 1 unmodified line
730	926	Math_y1(vm_registers_t &r, vm_t &vm)
731	927	{
732	928	auto &fp = r.stack.peek_frame()->base<F_Math_y1>();
	929	+
	930	+ // [PAL] Consume platform API
733	931	throw vm_system_exception{ "Function not implemented" };
734	932	}
735	933
@@ -737,5 +935,7 @@ void · 1 unmodified line
737	935	Math_yn(vm_registers_t &r, vm_t &vm)
738	936	{
739	937	auto &fp = r.stack.peek_frame()->base<F_Math_yn>();
	938	+
	939	+ // [PAL] Consume platform API
740	940	throw vm_system_exception{ "Function not implemented" };
741	941	}

src/vm/math/dgemm.c+411−0

@@ -0,0 +1,411 @@
	1	+
	2	+/* -- translated by f2c (version 19940927).
	3	+ You must link the resulting object file with the libraries:
	4	+ -lf2c -lm (in that order)
	5	+*/
	6	+
	7	+//#include "f2c.h"
	8	+
	9	+typedef int integer;
	10	+typedef int logical;
	11	+typedef double doublereal;
	12	+
	13	+#define max(a,b) ((a) > (b) ? (a) : (b))
	14	+
	15	+// Returns '1' based index of invalid argument, otherwise '0'
	16	+/* Subroutine / int dgemm_(char transa, char transb, integer m, integer *
	17	+ n, integer k, doublereal alpha, doublereal a, integer lda,
	18	+ doublereal b, integer ldb, doublereal beta, doublereal c, integer
	19	+ *ldc, int a_len, int b_len, int c_len)
	20	+{
	21	+
	22	+
	23	+ /* System generated locals */
	24	+ integer /a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, / i__1, i__2,
	25	+ i__3;
	26	+
	27	+ /* Local variables */
	28	+ integer info;
	29	+ logical nota, notb;
	30	+ doublereal temp;
	31	+ integer i, j, l, ncola, ncolb;
	32	+ extern logical lsame_(char , char );
	33	+ integer nrowa, nrowb;
	34	+ //extern /* Subroutine / int xerbla_(char , integer *);
	35	+
	36	+
	37	+/* Purpose
	38	+ =======
	39	+
	40	+ DGEMM performs one of the matrix-matrix operations
	41	+
	42	+ C := alphaop( A )op( B ) + beta*C,
	43	+
	44	+ where op( X ) is one of
	45	+
	46	+ op( X ) = X or op( X ) = X',
	47	+
	48	+ alpha and beta are scalars, and A, B and C are matrices, with op( A )
	49	+
	50	+ an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
	51	+
	52	+
	53	+ Parameters
	54	+ ==========
	55	+
	56	+ TRANSA - CHARACTER*1.
	57	+ On entry, TRANSA specifies the form of op( A ) to be used in
	58	+
	59	+ the matrix multiplication as follows:
	60	+
	61	+ TRANSA = 'N' or 'n', op( A ) = A.
	62	+
	63	+ TRANSA = 'T' or 't', op( A ) = A'.
	64	+
	65	+ TRANSA = 'C' or 'c', op( A ) = A'.
	66	+
	67	+ Unchanged on exit.
	68	+
	69	+ TRANSB - CHARACTER*1.
	70	+ On entry, TRANSB specifies the form of op( B ) to be used in
	71	+
	72	+ the matrix multiplication as follows:
	73	+
	74	+ TRANSB = 'N' or 'n', op( B ) = B.
	75	+
	76	+ TRANSB = 'T' or 't', op( B ) = B'.
	77	+
	78	+ TRANSB = 'C' or 'c', op( B ) = B'.
	79	+
	80	+ Unchanged on exit.
	81	+
	82	+ M - INTEGER.
	83	+ On entry, M specifies the number of rows of the matrix
	84	+
	85	+ op( A ) and of the matrix C. M must be at least zero.
	86	+
	87	+ Unchanged on exit.
	88	+
	89	+ N - INTEGER.
	90	+ On entry, N specifies the number of columns of the matrix
	91	+
	92	+ op( B ) and the number of columns of the matrix C. N must be
	93	+
	94	+ at least zero.
	95	+ Unchanged on exit.
	96	+
	97	+ K - INTEGER.
	98	+ On entry, K specifies the number of columns of the matrix
	99	+
	100	+ op( A ) and the number of rows of the matrix op( B ). K must
	101	+
	102	+ be at least zero.
	103	+ Unchanged on exit.
	104	+
	105	+ ALPHA - DOUBLE PRECISION.
	106	+ On entry, ALPHA specifies the scalar alpha.
	107	+ Unchanged on exit.
	108	+
	109	+ A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
	110	+
	111	+ k when TRANSA = 'N' or 'n', and is m otherwise.
	112	+ Before entry with TRANSA = 'N' or 'n', the leading m by k
	113	+
	114	+ part of the array A must contain the matrix A, otherwise
	115	+
	116	+ the leading k by m part of the array A must contain the
	117	+
	118	+ matrix A.
	119	+ Unchanged on exit.
	120	+
	121	+ LDA - INTEGER.
	122	+ On entry, LDA specifies the first dimension of A as declared
	123	+
	124	+ in the calling (sub) program. When TRANSA = 'N' or 'n' then
	125	+
	126	+ LDA must be at least max( 1, m ), otherwise LDA must be at
	127	+
	128	+ least max( 1, k ).
	129	+ Unchanged on exit.
	130	+
	131	+ B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
	132	+
	133	+ n when TRANSB = 'N' or 'n', and is k otherwise.
	134	+ Before entry with TRANSB = 'N' or 'n', the leading k by n
	135	+
	136	+ part of the array B must contain the matrix B, otherwise
	137	+
	138	+ the leading n by k part of the array B must contain the
	139	+
	140	+ matrix B.
	141	+ Unchanged on exit.
	142	+
	143	+ LDB - INTEGER.
	144	+ On entry, LDB specifies the first dimension of B as declared
	145	+
	146	+ in the calling (sub) program. When TRANSB = 'N' or 'n' then
	147	+
	148	+ LDB must be at least max( 1, k ), otherwise LDB must be at
	149	+
	150	+ least max( 1, n ).
	151	+ Unchanged on exit.
	152	+
	153	+ BETA - DOUBLE PRECISION.
	154	+ On entry, BETA specifies the scalar beta. When BETA is
	155	+
	156	+ supplied as zero then C need not be set on input.
	157	+ Unchanged on exit.
	158	+
	159	+ C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
	160	+ Before entry, the leading m by n part of the array C must
	161	+
	162	+ contain the matrix C, except when beta is zero, in which
	163	+
	164	+ case C need not be set on entry.
	165	+ On exit, the array C is overwritten by the m by n matrix
	166	+
	167	+ ( alphaop( A )op( B ) + beta*C ).
	168	+
	169	+ LDC - INTEGER.
	170	+ On entry, LDC specifies the first dimension of C as declared
	171	+
	172	+ in the calling (sub) program. LDC must be at least
	173	+
	174	+ max( 1, m ).
	175	+ Unchanged on exit.
	176	+
	177	+
	178	+ Level 3 Blas routine.
	179	+
	180	+ -- Written on 8-February-1989.
	181	+ Jack Dongarra, Argonne National Laboratory.
	182	+ Iain Duff, AERE Harwell.
	183	+ Jeremy Du Croz, Numerical Algorithms Group Ltd.
	184	+ Sven Hammarling, Numerical Algorithms Group Ltd.
	185	+
	186	+
	187	+
	188	+ Set NOTA and NOTB as true if A and B respectively are not
	189	+
	190	+ transposed and set NROWA, NCOLA and NROWB as the number of rows
	191	+
	192	+ and columns of A and the number of rows of B respectively.
	193	+
	194	+
	195	+
	196	+ Parameter adjustments
	197	+ Function Body */
	198	+
	199	+#define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]
	200	+#define B(I,J) b[(I)-1 + ((J)-1)* ( *ldb)]
	201	+#define C(I,J) c[(I)-1 + ((J)-1)* ( *ldc)]
	202	+
	203	+ nota = lsame_(transa, "N");
	204	+ notb = lsame_(transb, "N");
	205	+ if (nota) {
	206	+ nrowa = *m;
	207	+ ncola = *k;
	208	+ } else {
	209	+ nrowa = *k;
	210	+ ncola = *m;
	211	+ }
	212	+ if (notb) {
	213	+ nrowb = *k;
	214	+ ncolb = *n;
	215	+ } else {
	216	+ nrowb = *n;
	217	+ ncolb = *k;
	218	+ }
	219	+
	220	+/* Test the input parameters. */
	221	+
	222	+ info = 0;
	223	+ if (! nota && ! lsame_(transa, "C") && ! lsame_(transa, "T")) {
	224	+ info = 1;
	225	+ } else if (! notb && ! lsame_(transb, "C") && ! lsame_(transb,
	226	+ "T")) {
	227	+ info = 2;
	228	+ } else if (*m < 0) {
	229	+ info = 3;
	230	+ } else if (*n < 0) {
	231	+ info = 4;
	232	+ } else if (*k < 0) {
	233	+ info = 5;
	234	+ } else if (((lda) (ncola - 1)) > a_len) {
	235	+ info = 7; // A matrix
	236	+ } else if (*lda < max(1,nrowa)) {
	237	+ info = 8;
	238	+ } else if (((ldb) (ncolb - 1)) > b_len) {
	239	+ info = 9; // B matrix
	240	+ } else if (*ldb < max(1,nrowb)) {
	241	+ info = 10;
	242	+ } else if (((ldc) (*n - 1)) > c_len) {
	243	+ info = 12; // C matrix
	244	+ } else if (ldc < max(1,m)) {
	245	+ info = 13;
	246	+ }
	247	+
	248	+ if (info != 0) {
	249	+ //xerbla_("DGEMM ", &info);
	250	+ return info;
	251	+ }
	252	+
	253	+/* Quick return if possible. */
	254	+
	255	+ if (m == 0 \|\| n == 0 \|\| (alpha == 0. \|\| k == 0) && *beta == 1.) {
	256	+ return 0;
	257	+ }
	258	+
	259	+/* And if alpha.eq.zero. */
	260	+
	261	+ if (*alpha == 0.) {
	262	+ if (*beta == 0.) {
	263	+ i__1 = *n;
	264	+ for (j = 1; j <= *n; ++j) {
	265	+ i__2 = *m;
	266	+ for (i = 1; i <= *m; ++i) {
	267	+ C(i,j) = 0.;
	268	+/* L10: */
	269	+ }
	270	+/* L20: */
	271	+ }
	272	+ } else {
	273	+ i__1 = *n;
	274	+ for (j = 1; j <= *n; ++j) {
	275	+ i__2 = *m;
	276	+ for (i = 1; i <= *m; ++i) {
	277	+ C(i,j) = beta C(i,j);
	278	+/* L30: */
	279	+ }
	280	+/* L40: */
	281	+ }
	282	+ }
	283	+ return 0;
	284	+ }
	285	+
	286	+/* Start the operations. */
	287	+
	288	+ if (notb) {
	289	+ if (nota) {
	290	+
	291	+/* Form C := alphaAB + betaC. /
	292	+
	293	+ i__1 = *n;
	294	+ for (j = 1; j <= *n; ++j) {
	295	+ if (*beta == 0.) {
	296	+ i__2 = *m;
	297	+ for (i = 1; i <= *m; ++i) {
	298	+ C(i,j) = 0.;
	299	+/* L50: */
	300	+ }
	301	+ } else if (*beta != 1.) {
	302	+ i__2 = *m;
	303	+ for (i = 1; i <= *m; ++i) {
	304	+ C(i,j) = beta C(i,j);
	305	+/* L60: */
	306	+ }
	307	+ }
	308	+ i__2 = *k;
	309	+ for (l = 1; l <= *k; ++l) {
	310	+ if (B(l,j) != 0.) {
	311	+ temp = alpha B(l,j);
	312	+ i__3 = *m;
	313	+ for (i = 1; i <= *m; ++i) {
	314	+ C(i,j) += temp * A(i,l);
	315	+/* L70: */
	316	+ }
	317	+ }
	318	+/* L80: */
	319	+ }
	320	+/* L90: */
	321	+ }
	322	+ } else {
	323	+
	324	+/* Form C := alphaA'B + betaC /
	325	+
	326	+ i__1 = *n;
	327	+ for (j = 1; j <= *n; ++j) {
	328	+ i__2 = *m;
	329	+ for (i = 1; i <= *m; ++i) {
	330	+ temp = 0.;
	331	+ i__3 = *k;
	332	+ for (l = 1; l <= *k; ++l) {
	333	+ temp += A(l,i) * B(l,j);
	334	+/* L100: */
	335	+ }
	336	+ if (*beta == 0.) {
	337	+ C(i,j) = alpha temp;
	338	+ } else {
	339	+ C(i,j) = alpha temp + beta C(i,j);
	340	+ }
	341	+/* L110: */
	342	+ }
	343	+/* L120: */
	344	+ }
	345	+ }
	346	+ } else {
	347	+ if (nota) {
	348	+
	349	+/* Form C := alphaAB' + betaC /
	350	+
	351	+ i__1 = *n;
	352	+ for (j = 1; j <= *n; ++j) {
	353	+ if (*beta == 0.) {
	354	+ i__2 = *m;
	355	+ for (i = 1; i <= *m; ++i) {
	356	+ C(i,j) = 0.;
	357	+/* L130: */
	358	+ }
	359	+ } else if (*beta != 1.) {
	360	+ i__2 = *m;
	361	+ for (i = 1; i <= *m; ++i) {
	362	+ C(i,j) = beta C(i,j);
	363	+/* L140: */
	364	+ }
	365	+ }
	366	+ i__2 = *k;
	367	+ for (l = 1; l <= *k; ++l) {
	368	+ if (B(j,l) != 0.) {
	369	+ temp = alpha B(j,l);
	370	+ i__3 = *m;
	371	+ for (i = 1; i <= *m; ++i) {
	372	+ C(i,j) += temp * A(i,l);
	373	+/* L150: */
	374	+ }
	375	+ }
	376	+/* L160: */
	377	+ }
	378	+/* L170: */
	379	+ }
	380	+ } else {
	381	+
	382	+/* Form C := alphaA'B' + betaC /
	383	+
	384	+ i__1 = *n;
	385	+ for (j = 1; j <= *n; ++j) {
	386	+ i__2 = *m;
	387	+ for (i = 1; i <= *m; ++i) {
	388	+ temp = 0.;
	389	+ i__3 = *k;
	390	+ for (l = 1; l <= *k; ++l) {
	391	+ temp += A(l,i) * B(j,l);
	392	+/* L180: */
	393	+ }
	394	+ if (*beta == 0.) {
	395	+ C(i,j) = alpha temp;
	396	+ } else {
	397	+ C(i,j) = alpha temp + beta C(i,j);
	398	+ }
	399	+/* L190: */
	400	+ }
	401	+/* L200: */
	402	+ }
	403	+ }
	404	+ }
	405	+
	406	+ return 0;
	407	+
	408	+/* End of DGEMM . */
	409	+
	410	+} /* dgemm_ */
	411	+

src/vm/math/lsame.c+113−0

@@ -0,0 +1,113 @@
	1	+/* -- translated by f2c (version 19940927).
	2	+ You must link the resulting object file with the libraries:
	3	+ -lf2c -lm (in that order)
	4	+*/
	5	+
	6	+//#include "f2c.h"
	7	+
	8	+typedef int integer;
	9	+typedef int logical;
	10	+
	11	+logical lsame_(char ca, char cb)
	12	+{
	13	+
	14	+
	15	+ /* System generated locals */
	16	+ logical ret_val;
	17	+
	18	+ /* Local variables */
	19	+ integer inta, intb, zcode;
	20	+
	21	+
	22	+/* -- LAPACK auxiliary routine (version 2.0) --
	23	+ Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
	24	+ Courant Institute, Argonne National Lab, and Rice University
	25	+ January 31, 1994
	26	+
	27	+
	28	+ Purpose
	29	+ =======
	30	+
	31	+ LSAME returns .TRUE. if CA is the same letter as CB regardless of
	32	+ case.
	33	+
	34	+ Arguments
	35	+ =========
	36	+
	37	+ CA (input) CHARACTER*1
	38	+ CB (input) CHARACTER*1
	39	+ CA and CB specify the single characters to be compared.
	40	+
	41	+ =====================================================================
	42	+
	43	+
	44	+
	45	+ Test if the characters are equal */
	46	+
	47	+ ret_val = (unsigned char )ca == (unsigned char )cb;
	48	+ if (ret_val) {
	49	+ return ret_val;
	50	+ }
	51	+
	52	+/* Now test for equivalence if both characters are alphabetic. */
	53	+
	54	+ zcode = 'Z';
	55	+
	56	+/* Use 'Z' rather than 'A' so that ASCII can be detected on Prime
	57	+ machines, on which ICHAR returns a value with bit 8 set.
	58	+ ICHAR('A') on Prime machines returns 193 which is the same as
	59	+ ICHAR('A') on an EBCDIC machine. */
	60	+
	61	+ inta = (unsigned char )ca;
	62	+ intb = (unsigned char )cb;
	63	+
	64	+ if (zcode == 90 \|\| zcode == 122) {
	65	+
	66	+/* ASCII is assumed - ZCODE is the ASCII code of either lower o
	67	+r
	68	+ upper case 'Z'. */
	69	+
	70	+ if (inta >= 97 && inta <= 122) {
	71	+ inta += -32;
	72	+ }
	73	+ if (intb >= 97 && intb <= 122) {
	74	+ intb += -32;
	75	+ }
	76	+
	77	+ } else if (zcode == 233 \|\| zcode == 169) {
	78	+
	79	+/* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower
	80	+ or
	81	+ upper case 'Z'. */
	82	+
	83	+ if (inta >= 129 && inta <= 137 \|\| inta >= 145 && inta <= 153 \|\| inta
	84	+ >= 162 && inta <= 169) {
	85	+ inta += 64;
	86	+ }
	87	+ if (intb >= 129 && intb <= 137 \|\| intb >= 145 && intb <= 153 \|\| intb
	88	+ >= 162 && intb <= 169) {
	89	+ intb += 64;
	90	+ }
	91	+
	92	+ } else if (zcode == 218 \|\| zcode == 250) {
	93	+
	94	+/* ASCII is assumed, on Prime machines - ZCODE is the ASCII cod
	95	+e
	96	+ plus 128 of either lower or upper case 'Z'. */
	97	+
	98	+ if (inta >= 225 && inta <= 250) {
	99	+ inta += -32;
	100	+ }
	101	+ if (intb >= 225 && intb <= 250) {
	102	+ intb += -32;
	103	+ }
	104	+ }
	105	+ ret_val = inta == intb;
	106	+
	107	+/* RETURN
	108	+
	109	+ End of LSAME */
	110	+
	111	+ return ret_val;
	112	+} /* lsame_ */
	113	+