asps.c 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446
  1. /* file asps.c */
  2. /*
  3. * This code may be used and modified, and redistributed in binary
  4. * or source form, subject to the "CCL Public License", which should
  5. * accompany it. This license is a variant on the BSD license, and thus
  6. * permits use of code derived from this in either open and commercial
  7. * projects: but it does require that updates to this code be made
  8. * available back to the originators of the package.
  9. * Before merging other code in with this or linking this code
  10. * with other packages or libraries please check that the license terms
  11. * of the other material are compatible with those of this.
  12. */
  13. /* Signature: 6a193f20 08-Apr-2002 */
  14. #include <stdarg.h>
  15. #include <sys/types.h>
  16. #include "machine.h"
  17. #include "tags.h"
  18. #include "cslerror.h"
  19. #include "externs.h"
  20. #include "entries.h"
  21. #include "arith.h"
  22. #ifndef WINDOWS_NT
  23. #define __stdcall
  24. #endif
  25. /*
  26. * The overall design for the asps mechanism relies on the following facts:
  27. * 1. No NAG routine uses the same type of ASP more than once.
  28. * 2. We have a hard-coded limit of 10 Asps in use at any one time (in
  29. * theory we could make this 10 Asps per NAG routine by changing the
  30. * storage strategy for the code/environment pairs).
  31. * 3. If NAG routines call one another, they do it via the top-level AXIOM
  32. * interfaces (and so go through the setup/restore process.
  33. *
  34. * The code and environment are the values obtained by applying car and cdr
  35. * respectively to a spadclosure. They are used in calling funcall in lisp
  36. * as follows:
  37. * (FUNCALL code arg1 arg2 ... argn environment)
  38. *
  39. * When we call setup on an asp, we push the code part onto user_base_0 and
  40. * the environment onto user_base_1 (which we treat as Lisp lists), and then
  41. * store the index in user_base_0/user_base_1 into the appropriate variable
  42. * (this is the index counting from the base).
  43. *
  44. * When we call an asp we copy the appropriate code/environment portions onto
  45. * users_base_2 and user_base_3, and pass them as arguments to Lfuncall and
  46. * friends.
  47. *
  48. * When we call restore on an asp we pop user_base_0/user_base_1 (this is
  49. * dodgey in one sense, but we assume that we restore all the asps from a
  50. * NAG call at the same time so the order in which we pop the user_bases is
  51. * in fact immaterial). Finally we restore the old asp_n_pointer value.
  52. *
  53. */
  54. extern Lisp_Object user_base_0, user_base_1, user_base_2, user_base_3;
  55. #define CODE_STACK user_base_0
  56. #define ENVIRONMENT_STACK user_base_1
  57. #define CODE user_base_2
  58. #define ENVIRONMENT user_base_3
  59. /* The number of asps on the stacks at present. */
  60. int32 stack_depth=0;
  61. /*
  62. * The current location of an asp's code/environment is stored on a stack, so
  63. * that we can nest calls which use asps.
  64. *
  65. */
  66. typedef struct Asp_Loc {
  67. int32 loc;
  68. struct Asp_Loc *next;
  69. } Asp_Loc;
  70. Asp_Loc *asp1_index=NULL,*asp4_index=NULL,*asp6_index=NULL,*asp35_index=NULL;
  71. #define location(asp) \
  72. (asp == NULL) ? aerror0("Asp calling mechanism corrupted"): asp -> loc;
  73. Asp_Loc *push_asp_num(Asp_Loc *asp) {
  74. Asp_Loc *new;
  75. new = (Asp_Loc *)malloc(sizeof(Asp_Loc));
  76. new->next = asp;
  77. new->loc = stack_depth;
  78. return(new);
  79. }
  80. #define pop_asp_num(asp) asp->next
  81. void push_asp(Lisp_Object nil, Lisp_Object f_code, Lisp_Object f_env)
  82. {
  83. Lisp_Object w;
  84. if (stack_depth == 0) {
  85. push(f_env);
  86. w = ncons(f_code);
  87. pop(f_env);
  88. errexitv();
  89. f_code = w;
  90. push(f_code);
  91. w = ncons(f_env);
  92. pop(f_code);
  93. errexitv();
  94. f_env = w;
  95. CODE_STACK = f_code;
  96. ENVIRONMENT_STACK = f_env;
  97. }
  98. else {
  99. push(f_env);
  100. w = Lcons(nil,f_code,CODE_STACK);
  101. pop(f_env);
  102. errexitv();
  103. f_code = w;
  104. push(f_code);
  105. w = Lcons(nil,f_env,ENVIRONMENT_STACK);
  106. pop(f_code);
  107. errexitv();
  108. f_env = w;
  109. CODE_STACK = f_code;
  110. ENVIRONMENT_STACK = f_env;
  111. }
  112. ++ stack_depth;
  113. }
  114. void pop_asp(Lisp_Object nil)
  115. {
  116. CODE_STACK=qcdr(CODE_STACK);
  117. ENVIRONMENT_STACK=qcdr(ENVIRONMENT_STACK);
  118. -- stack_depth;
  119. }
  120. #define get_code(index) Lelt(C_nil,CODE_STACK,fixnum_of_int(stack_depth-index))
  121. #define get_env(index) Lelt(C_nil,ENVIRONMENT_STACK,fixnum_of_int(stack_depth-index))
  122. /*
  123. * This function restores the global variables which are used in calling asps
  124. * to their previous state.
  125. */
  126. Lisp_Object asp_restore(Lisp_Object nil, int32 asp)
  127. {
  128. pop_asp(nil);
  129. switch (int_of_fixnum(asp))
  130. {
  131. case 1:
  132. asp1_index = pop_asp_num(asp1_index);
  133. break;
  134. case 4:
  135. asp4_index = pop_asp_num(asp4_index);
  136. break;
  137. case 6:
  138. asp6_index = pop_asp_num(asp6_index);
  139. break;
  140. case 35:
  141. asp35_index = pop_asp_num(asp35_index);
  142. break;
  143. default:
  144. aerror1("asp_restore: Unknown asp type:",asp);
  145. break;
  146. }
  147. return onevalue(lisp_true);
  148. }
  149. /*
  150. * This function sets up the global variables which are needed in calling asps
  151. */
  152. Lisp_Object MS_CDECL asp_set(Lisp_Object nil, int nargs, ...)
  153. {
  154. int32 asp;
  155. va_list args;
  156. Lisp_Object aspnum, f_car, f_cdr;
  157. argcheck(nargs,3,"asp_set");
  158. va_start(args,nargs);
  159. aspnum = va_arg(args, Lisp_Object);
  160. f_car = va_arg(args, Lisp_Object);
  161. f_cdr = va_arg(args, Lisp_Object);
  162. va_end(args);
  163. asp = int_of_fixnum(aspnum);
  164. switch (asp)
  165. {
  166. case 1:
  167. push_asp(nil,f_car, f_cdr);
  168. errexit();
  169. asp1_index = push_asp_num(asp1_index);
  170. break;
  171. case 4:
  172. push_asp(nil,f_car, f_cdr);
  173. errexit();
  174. asp4_index = push_asp_num(asp4_index);
  175. break;
  176. case 6:
  177. push_asp(nil,f_car, f_cdr);
  178. errexit();
  179. asp6_index = push_asp_num(asp6_index);
  180. break;
  181. case 35:
  182. push_asp(nil,f_car, f_cdr);
  183. errexit();
  184. asp35_index = push_asp_num(asp35_index);
  185. break;
  186. default:
  187. aerror1("asp_set: Unknown asp type:",aspnum);
  188. break;
  189. }
  190. return onevalue(lisp_true);
  191. }
  192. /* Routines for converting between representations. */
  193. Lisp_Object mkAXIOMVectorDF(double *v, int32 dim)
  194. {
  195. Lisp_Object nil=C_nil;
  196. Lisp_Object new, Lflt;
  197. int32 i;
  198. new = getvector(TAG_VECTOR, TYPE_SIMPLE_VEC, 4*dim+4);
  199. errexit();
  200. /* vectors must pad to an even number of words */
  201. if ((dim & 1) == 0) elt(new,dim) = nil;
  202. for (i=0;i<dim;++i) {
  203. push(new);
  204. Lflt = make_boxfloat(*(v+i),TYPE_DOUBLE_FLOAT);
  205. pop(new);
  206. errexit();
  207. elt(new,i) = Lflt;
  208. }
  209. return onevalue(new);
  210. }
  211. void mkFortranVectorDouble(double *loc, Lisp_Object v, int32 dim)
  212. {
  213. int32 i;
  214. /*Lisp_Object nil=C_nil;*/
  215. for (i=0;i<dim;++i) {
  216. push(v);
  217. *(loc + i) = float_of_number(elt(v,i));
  218. pop(v);
  219. }
  220. }
  221. void mkFortran2DArrayDouble(double *loc, Lisp_Object v, int32 rows, int32 cols)
  222. {
  223. int32 i,j;
  224. /*Lisp_Object nil=C_nil;*/
  225. for (i=0;i<rows;++i)
  226. for (j=0;j<cols;++j) {
  227. push(v);
  228. *(loc + j*cols + i) = float_of_number(elt(v,i*cols+j));
  229. pop(v);
  230. }
  231. }
  232. /* Code for ASP1 */
  233. double __stdcall asp1 (double *x)
  234. {
  235. Lisp_Object arg = make_boxfloat(*x,TYPE_DOUBLE_FLOAT), result, code, env;
  236. Lisp_Object nil = C_nil;
  237. code = get_code(asp1_index->loc);
  238. env = get_env(asp1_index->loc);
  239. result = Lfuncalln(C_nil,3,code,arg,env);
  240. errexit();
  241. if (exception_pending()) aerror0("Error in evaluating function (ASP1)");
  242. return float_of_number(result);
  243. }
  244. #ifdef TEST_ASPS
  245. Lisp_Object test_asp1(Lisp_Object nil, Lisp_Object a)
  246. {
  247. double arg = float_of_number(a);
  248. arg=asp1(&arg);
  249. return make_boxfloat(arg,TYPE_DOUBLE_FLOAT);
  250. }
  251. #endif
  252. /* Code for ASP4 */
  253. double __stdcall asp4 (int32 *ndim, double *x)
  254. {
  255. Lisp_Object arg, result, code, env;
  256. arg = mkAXIOMVectorDF(x,*ndim);
  257. code = get_code(asp4_index->loc);
  258. env = get_env(asp4_index->loc);
  259. result = Lfuncalln(C_nil,3,code,arg,env);
  260. return float_of_number(result);
  261. }
  262. #ifdef TEST_ASPS
  263. Lisp_Object test_asp4(Lisp_Object nil, Lisp_Object v)
  264. {
  265. double *x;
  266. int ndim;
  267. Lisp_Object result;
  268. ndim = (length_of_header(vechdr(v)) - 4)/4;
  269. x = (double *)malloc(ndim*sizeof(double));
  270. mkFortranVectorDouble(x,v,ndim);
  271. result = make_boxfloat(asp4(&ndim,x),TYPE_DOUBLE_FLOAT);
  272. free(x);
  273. return(result);
  274. }
  275. #endif
  276. /* Code for ASP6 */
  277. void __stdcall asp6 (int32 *n, double *x, double *fvec, int32 *iflag)
  278. {
  279. Lisp_Object arg, code, env;
  280. arg = mkAXIOMVectorDF(x,*n);
  281. code = get_code(asp6_index->loc);
  282. env = get_env(asp6_index->loc);
  283. mkFortranVectorDouble(fvec,Lfuncalln(C_nil,3,code,arg,env),*n);
  284. }
  285. #ifdef TEST_ASPS
  286. Lisp_Object test_asp6(Lisp_Object nil, Lisp_Object v)
  287. {
  288. double *x, *fvec;
  289. int n,iflag;
  290. Lisp_Object result;
  291. n = (length_of_header(vechdr(v)) - 4)/4;
  292. x = (double *)malloc(n*sizeof(double));
  293. fvec = (double *)malloc(n*sizeof(double));
  294. mkFortranVectorDouble(x,v,n);
  295. asp6(&n,x,fvec,&iflag);
  296. result = mkAXIOMVectorDF(fvec,n);
  297. free(x);
  298. free(fvec);
  299. return(result);
  300. }
  301. #endif
  302. /* Code for ASP35 */
  303. void __stdcall asp35 (int32 *n, double *x, double *fvec, double *fjac,
  304. int32 *ldfjac, int32 *iflag)
  305. {
  306. Lisp_Object Lx, Liflag, Lresult, code, env;
  307. Lx = mkAXIOMVectorDF(x,*n);
  308. Liflag = fixnum_of_int(*iflag);
  309. code = get_code(asp35_index->loc);
  310. env = get_env(asp35_index->loc);
  311. Lresult = Lfuncalln(C_nil,4,code,Lx,Liflag,env);
  312. if (*iflag == 1)
  313. mkFortranVectorDouble(fvec,Lresult,*n);
  314. else
  315. mkFortran2DArrayDouble(fjac,Lresult,*ldfjac,*n);
  316. }
  317. #ifdef TEST_ASPS
  318. Lisp_Object test_asp35(Lisp_Object nil, Lisp_Object v, Lisp_Object flag)
  319. {
  320. double *x, *fvec, *fjac;
  321. int n, ldfjac, iflag;
  322. Lisp_Object result;
  323. n = (length_of_header(vechdr(v)) - 4)/4;
  324. ldfjac=n;
  325. x = (double *)malloc(n*sizeof(double));
  326. mkFortranVectorDouble(x,v,n);
  327. fjac = (double *)malloc(n*ldfjac*sizeof(double));
  328. fvec = (double *)malloc(n*sizeof(double));
  329. iflag=int_of_fixnum(flag);
  330. asp35(&n,x,fvec,fjac,&ldfjac,&iflag);
  331. if (iflag == 1)
  332. result = mkAXIOMVectorDF(fvec,n);
  333. else
  334. result = mkAXIOMVectorDF(fjac,n*ldfjac);
  335. free(x);
  336. free(fjac);
  337. free(fvec);
  338. return(result);
  339. }
  340. #endif
  341. #ifndef TEST_ASPS
  342. Lisp_Object asp_error1(Lisp_Object nil, Lisp_Object a1)
  343. {
  344. return aerror0("The Windows version of the NAG Link is not installed");
  345. }
  346. Lisp_Object asp_error2(Lisp_Object nil, Lisp_Object a1, Lisp_Object a2)
  347. {
  348. return aerror0("The Windows version of the NAG Link is not installed");
  349. }
  350. Lisp_Object MS_CDECL asp_error0(Lisp_Object nil, int32 nargs, ...)
  351. {
  352. return aerror0("The Windows version of the NAG Link is not installed");
  353. }
  354. #endif
  355. setup_type const asp_setup[] =
  356. {
  357. {"asp-setup", wrong_no_na, wrong_no_nb, asp_set},
  358. {"asp-restore", asp_restore, too_many_1, wrong_no_1},
  359. #ifdef TEST_ASPS
  360. {"asp1", test_asp1, too_many_1, wrong_no_1},
  361. {"asp4", test_asp4, too_many_1, wrong_no_1},
  362. {"asp6", test_asp6, too_many_1, wrong_no_1},
  363. {"asp35", too_few_2, test_asp35, wrong_no_2},
  364. #else
  365. {"asp1", asp_error1, asp_error2, asp_error0},
  366. {"asp4", asp_error1, asp_error2, asp_error0},
  367. {"asp6", asp_error1, asp_error2, asp_error0},
  368. {"asp35", asp_error1, asp_error2, asp_error0},
  369. #endif
  370. {NULL, 0, 0, 0}
  371. };
  372. /* end of file asps.c */