1 /*
2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
3 * Applied Mathematics, Norway.
4 *
5 * Contact information: E-mail: tor.dokken@sintef.no
6 * SINTEF ICT, Department of Applied Mathematics,
7 * P.O. Box 124 Blindern,
8 * 0314 Oslo, Norway.
9 *
10 * This file is part of SISL.
11 *
12 * SISL is free software: you can redistribute it and/or modify
13 * it under the terms of the GNU Affero General Public License as
14 * published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * SISL is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Affero General Public License for more details.
21 *
22 * You should have received a copy of the GNU Affero General Public
23 * License along with SISL. If not, see
24 * <http://www.gnu.org/licenses/>.
25 *
26 * In accordance with Section 7(b) of the GNU Affero General Public
27 * License, a covered work must retain the producer line in every data
28 * file that is created or manipulated using SISL.
29 *
30 * Other Usage
31 * You can be released from the requirements of the license by purchasing
32 * a commercial license. Buying such a license is mandatory as soon as you
33 * develop commercial activities involving the SISL library without
34 * disclosing the source code of your own applications.
35 *
36 * This file may be used in accordance with the terms contained in a
37 * written agreement between you and SINTEF ICT.
38 */
39
40 #include "sisl-copyright.h"
41
42 /*
43 *
44 * $Id: s1700.c,v 1.1 1994-04-21 12:10:42 boh Exp $
45 *
46 */
47
48
49 #define S1700
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s1700(int imy,int ik,int in,int iv,int * jpl,int * jfi,int * jla,double * et,double apar,double * galfa,int * jstat)55 s1700(int imy,int ik,int in,int iv,
56 int *jpl,int *jfi,int *jla,double *et,double apar,double *galfa,int *jstat)
57 #else
58 void s1700(imy,ik,in,iv,jpl,jfi,jla,et,apar,galfa,jstat)
59 int imy;
60 int ik;
61 int in;
62 int iv;
63 int *jpl;
64 int *jfi;
65 int *jla;
66 double *et;
67 double apar;
68 double *galfa;
69 int *jstat;
70 #endif
71 /*
72 ********************************************************************
73 *
74 *********************************************************************
75 *
76 * PURPOSE : To compute in a compact format a line in the discrete
77 * B-spline matrix converting between an orginal basis
78 * "et" and a new basis consisting of "et" and one p-tuppel
79 * knot, where 0 < p <= ik.
80 * To insert one p-tuppel knot the function must be called
81 * ik - (multiplicity of the knot) + 2*(p-1) times.
82 * The function can also be used for inserting knots at
83 * more than one value if at least ik-1 knots exist in "et"
84 * between the values in question.
85 * In this case apar must be changed and the index of the
86 * new vertices must be adjusted in the proper manner
87 * by the calling function.
88 *
89 *
90 *
91 * INPUT : imy - If one p-tuppel knot are to be inserted, and
92 * the index of the new vertice are j then
93 * imy = j when et[j] < apar and then unchanged
94 * for the rest of j.
95 * If knots are to be inserted in more than one
96 * value, index the new vertices as it is only
97 * one value to be insert, and adjust the index
98 * after calling this function.
99 * ik - The order of the B-spline.
100 * in - The number of the orginal vertices.
101 * iv - The number of new knots to insert in the
102 * area between knot number j and knot number j+ik
103 * in the new knot vector.
104 * et - The knot vector.
105 * apar - The value where to insert new knots.
106 *
107 *
108 *
109 * OUTPUT : jpl - The negativ difference between the index in galfa
110 * and the real knot inserten matrix.
111 * jfi - The index of the first element in the line j in the
112 * the real knot inserten matrix whice is not zero.
113 * The element with the index (jfi+jpl) in galfa
114 * is the same as the element with index jfi in
115 * the real line j in the knot inserten matrix.
116 * jla - The index of the last element in the line j in the
117 * real knot inserten matrix whice is not zero.
118 * The element with the index (jla+jpl) in galfa
119 * is the same as the element with index jla in
120 * the real line j in the knot inserten matrix.
121 * galfa - A compressed line in the knot inserten matrix.
122 * jstat - status messages
123 * > 0 : warning
124 * = 0 : ok
125 * < 0 : error
126 *
127 *
128 * METHOD : Using the Oslo-algorithm
129 *
130 *
131 * REFERENCES : Making The Oslo algorithm more efficient.
132 * by T.Lyche and K.Moerken.
133 * SIAM J.NUMER.ANAL Vol. 23, No. 3, June 1986.
134 *
135 *-
136 * CALLS :
137 *
138 * WRITTEN BY : Arne Laksaa, SI, 88-06.
139 *
140 **********************************************************************/
141 {
142 int kpos=0; /* Posisjon of error. */
143 int kj,kv; /* Help variable */
144 int kp; /* Control variable in loop. */
145 double *salfa; /* Help pointer to galfa. */
146 double tbeta,tbeta1; /* Help variabels */
147 double td1,td2; /* Help variabels */
148 double *t1,*t2; /* Pointers to the knot vector. */
149
150
151 /* Check that the number of knots we insert is not to large */
152
153 if (iv >= ik) goto err152;
154
155
156 /* Compute the negativ difference between the index in galfa and
157 the real knot inserten matrix. */
158
159 *jpl=ik-imy-1;
160
161
162 /* Changing the galfa so we may use the index in the real matrix. */
163
164 galfa += *jpl;
165
166
167 /* Initialise the last element. */
168
169 galfa[imy] = 1;
170
171
172 /* Here we go one time for each new knot we insert. */
173
174 for (kj=in+iv-2,in+=ik-1,kv=ik-iv,kp=0; kp<iv; kp++,kv++)
175 {
176 /* The initialising: The two first are not changing.
177 kj = in+iv-2, minus the maximum of kp it
178 gives the index of the last
179 orginal vertices.
180 in = in+ik-1, the index of the last element in et.
181 kv = ik-iv , the nuber of old knots in the field.
182 This variabel is counting up to ik
183 (the order) during the loops. */
184
185
186 /* Here we note the special case where we are at the
187 start of the matrix and we does not have a k-touple
188 knot at this end. */
189
190 if (kp>=imy) tbeta1=(apar - *et)* *galfa/(et[kv] - *et);
191 else tbeta1=(double)0.0;
192
193
194 *jfi=max(1,imy-kp); *jla=min(imy,kj-kp);
195
196
197 /* For details about this loop look in the reference. */
198
199 for (salfa=galfa+*jfi,t1=et+*jfi,t2=et+*jla; t1<=t2; t1++,salfa++)
200 {
201 td1 = apar - *t1;
202 td2 = t1[kv] - apar;
203 tbeta = *salfa/(td1 + td2);
204 salfa[-1] = td2*tbeta + tbeta1;
205 tbeta1 = td1*tbeta;
206 }
207
208
209 /* Here we note the special case where we are at the
210 end of the matrix and we does not have a k-touple
211 knot at this end. */
212
213 if (*jla<imy)
214 {
215 t1 = et + in;
216 *(salfa-1) = tbeta1+(*t1-apar)* *salfa/(*t1 - *(t2+1));
217 } else *(salfa-1) = tbeta1;
218 }
219
220
221 /* Adjusting the index of first and last in galfa. */
222
223 if (iv) (*jfi)--;
224 else *jfi = *jla = imy;
225
226
227 /* Updating output. */
228
229 *jstat = 0;
230 goto out;
231
232
233 /* Error, to many insertions knots. */
234
235 err152:
236 *jstat = -152;
237 s6err("s1700",*jstat,kpos);
238 goto out;
239
240 out:
241 return;
242 }
243