1 /* ======================================================================= */
2 /* DSPF_dp_lud_d.c -- lower upper decomposition using double precision */
3 /* */
4 /* */
5 /* Rev 1.0.0 */
6 /* */
7 /* Copyright (C) 2013 Texas Instruments Incorporated - http://www.ti.com/ */
8 /* */
9 /* */
10 /* Redistribution and use in source and binary forms, with or without */
11 /* modification, are permitted provided that the following conditions */
12 /* are met: */
13 /* */
14 /* Redistributions of source code must retain the above copyright */
15 /* notice, this list of conditions and the following disclaimer. */
16 /* */
17 /* Redistributions in binary form must rep roduce the above copyright */
18 /* notice, this list of conditions and the following disclaimer in the */
19 /* documentation and/or other materials provided with the */
20 /* distribution. */
21 /* */
22 /* Neither the name of Texas Instruments Incorporated nor the names of */
23 /* its contributors may be used to endorse or promote products derived */
24 /* from this software without specific prior written permission. */
25 /* */
26 /* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */
27 /* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */
28 /* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR */
29 /* A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT */
30 /* OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, */
31 /* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT */
32 /* LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, */
33 /* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY */
34 /* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT */
35 /* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE */
36 /* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
37 /* */
38 /* ======================================================================= */
39 #include <stdlib.h>
40 #include <stdio.h>
41 #include <math.h>
42 #include <c6x.h>
44 /* ======================================================================= */
45 /* Interface header files for the natural C and optimized C code */
46 /* ======================================================================= */
47 #include "DSPF_dp_lud_cn.h"
48 #include "DSPF_dp_lud.h"
50 /* ======================================================================= */
51 /* define size constants */
52 /* ======================================================================= */
53 #if defined(__TI_EABI__)
54 #define kernel_size _kernel_size
55 #endif
57 /* ======================================================================= */
58 /* extern parameters */
59 /* ======================================================================= */
60 extern char kernel_size;
62 /* ======================================================================= */
63 /* define constants */
64 /* ======================================================================= */
65 #define MAX_MATRIX_ORDER 128 /* 164 is max for L2SRAM */
66 #define MAX_MATRIX_SIZE (MAX_MATRIX_ORDER*MAX_MATRIX_ORDER)
68 /* ======================================================================= */
69 /* formula arrays */
70 /* ======================================================================= */
71 #define CYCLE_FORMULA_ORDER_PT1 64
72 #define CYCLE_FORMULA_ORDER_PT2 128
73 #define FORMULA_SIZE 2
74 long long form_cycle[FORMULA_SIZE];
76 /* ======================================================================= */
77 /* double word alignment of arrays */
78 /* ======================================================================= */
79 #pragma DATA_SECTION(A,".far:matrix")
80 #pragma DATA_ALIGN(A,8)
81 double A[MAX_MATRIX_SIZE];
83 #pragma DATA_SECTION(L,".far:matrix")
84 #pragma DATA_ALIGN(L,8)
85 double L[MAX_MATRIX_SIZE];
87 #pragma DATA_SECTION(U,".far:matrix")
88 #pragma DATA_ALIGN(U,8)
89 double U[MAX_MATRIX_SIZE];
91 #pragma DATA_SECTION(P,".far:matrix")
92 #pragma DATA_ALIGN(P,1)
93 unsigned short P[MAX_MATRIX_SIZE];
95 #pragma DATA_SECTION(LU,".far:matrix")
96 #pragma DATA_ALIGN(LU,8)
97 double LU[MAX_MATRIX_SIZE];
99 #pragma DATA_SECTION(inv_A,".far:matrix")
100 #pragma DATA_ALIGN(inv_A,8)
101 double inv_A[MAX_MATRIX_SIZE];
103 #pragma DATA_SECTION(b,".far:matrix")
104 #pragma DATA_ALIGN(b,8)
105 double b[MAX_MATRIX_ORDER];
107 #pragma DATA_SECTION(b_mod,".far:matrix")
108 #pragma DATA_ALIGN(b_mod,8)
109 double b_mod[MAX_MATRIX_ORDER];
111 #pragma DATA_SECTION(x,".far:matrix")
112 #pragma DATA_ALIGN(x,8)
113 double x[MAX_MATRIX_ORDER];
115 #pragma DATA_SECTION(y,".far:matrix")
116 #pragma DATA_ALIGN(y,8)
117 double y[MAX_MATRIX_ORDER];
119 /* ======================================================================= */
120 /* Main -- Top level driver for testing the algorithm */
121 /* ======================================================================= */
122 void main(void) {
124 unsigned short test,order,pass;
125 double sum,error;
126 int row,col,k,status;
127 long long t_overhead, t_start, t_stop;
128 long long lud_t_opt,lud_t_cn;
129 double tolerance=0.000001;
131 /* --------------------------------------------------------------------- */
132 /* Compute the overhead of calling clock twice to get timing info */
133 /* --------------------------------------------------------------------- */
134 TSCL=0, TSCH=0;
135 t_start = _itoll(TSCH,TSCL);
136 t_stop = _itoll(TSCH,TSCL);
137 t_overhead = t_stop - t_start;
139 /* --------------------------------------------------------------------- */
140 /* process test cases */
141 /* --------------------------------------------------------------------- */
142 pass=1;
143 for (test=1;test<=4;test++) {
145 printf("DSPF_dp_lud\tIter#: %d\t", test);
147 switch(test) {
148 case 1: { order=3;
149 A[0]=1;A[1]=2;A[2]=3;A[3]=4;A[4]=5;A[5]=6;A[6]=7;A[7]=8;A[8]=0; // non symmetric matrix
150 b[0]=1;b[1]=2;b[2]=3;
151 break;
152 }
153 case 2: { order=3;
154 A[0]=0;A[1]=1;A[2]=2;A[3]=3;A[4]=4;A[5]=5;A[6]=6;A[7]=7;A[8]=8; // non symmetric singular matrix
155 break;
156 }
157 case 3: { order=CYCLE_FORMULA_ORDER_PT1;
158 srand(1);
159 for (row=0;row<order;row++) b[row]=(double)row+1.0;
160 for (row=0;row<order;row++) { // random matrix
161 for (col=0;col<order;col++) {
162 A[row*order+col]=(double)(rand())/((double)RAND_MAX);
163 }
164 }
165 break;
166 }
167 case 4: { order=CYCLE_FORMULA_ORDER_PT2;
168 srand(1);
169 for (row=0;row<order;row++) b[row]=(double)row+1.0;
170 for (row=0;row<order;row++) { // random matrix
171 for (col=0;col<order;col++) {
172 A[row*order+col]=(double)(rand())/((double)RAND_MAX);
173 }
174 }
175 break;
176 }
177 }
179 /* ------------------------------------------------------------------- */
180 /* decompose A into P, L, and U where A=transpose(P)*L*U */
181 /* ------------------------------------------------------------------- */
182 TSCL=0, TSCH=0;
183 t_start = _itoll(TSCH,TSCL);
184 status = DSPF_dp_lud_cn(order,A,L,U,P);
185 t_stop = _itoll(TSCH,TSCL);
186 lud_t_cn = t_stop-t_start-t_overhead;
187 if (status==-1) {
188 printf("LU decomposition failed!\n");
189 }
191 /* ------------------------------------------------------------------- */
192 /* check decomposition: transpose(P)*L*U=A */
193 /* ------------------------------------------------------------------- */
194 for (row=0;row<order;row++) {
195 for (col=0;col<order;col++) {
196 sum=0;
197 for (k=0;k<order;k++) {
198 sum+=L[k+row*order]*U[col+k*order];
199 }
200 LU[col+row*order]=sum;
201 }
202 }
203 for (row=0;row<order;row++) {
204 for (col=0;col<order;col++) {
205 sum=0;
206 for (k=0;k<order;k++) {
207 sum+=P[row+k*order]*LU[col+k*order];
208 }
209 error=fabs(A[col+row*order]-sum);
210 if (error>tolerance) {
211 pass=0;
212 printf("error=%e\n",error);
213 }
214 }
215 }
217 /* ------------------------------------------------------------------- */
218 /* decompose A and generate L and U using optimized C code */
219 /* ------------------------------------------------------------------- */
220 TSCL=0, TSCH=0;
221 t_start = _itoll(TSCH,TSCL);
222 status=DSPF_dp_lud(order,A,L,U,P);
223 t_stop = _itoll(TSCH,TSCL);
224 if (status==-1) {
225 printf("LU decomposition failed!\n");
226 }
227 lud_t_opt=t_stop-t_start-t_overhead;
229 /* ------------------------------------------------------------------- */
230 /* check decomposition: transpose(P)*L*U=A */
231 /* ------------------------------------------------------------------- */
232 for (row=0;row<order;row++) {
233 for (col=0;col<order;col++) {
234 sum=0;
235 for (k=0;k<order;k++) {
236 sum+=L[k+row*order]*U[col+k*order];
237 }
238 LU[col+row*order]=sum;
239 }
240 }
241 for (row=0;row<order;row++) {
242 for (col=0;col<order;col++) {
243 sum=0;
244 for (k=0;k<order;k++) {
245 sum+=P[row+k*order]*LU[col+k*order];
246 }
247 error=fabs(A[col+row*order]-sum);
248 if (error>tolerance) {
249 pass=0;
250 printf("error=%e\n",error);
251 }
252 }
253 }
255 if (pass) {
256 printf("LUD Result Successful");
257 printf("\torder=%3d\tnatC: %lld\toptC: %lld\n",order, lud_t_cn, lud_t_opt);
258 } else {
259 printf("Result Failure\torder=%3d\n",order);
260 }
262 if (order==CYCLE_FORMULA_ORDER_PT1) {
263 form_cycle[0] = lud_t_opt;
264 } else if (order==CYCLE_FORMULA_ORDER_PT2) {
265 form_cycle[1] = lud_t_opt;
266 }
268 } /* for (test=1;test<=4;test++) */
270 /* ------------------------------------------------------------------- */
271 /* provide memory and cycles information */
272 /* ------------------------------------------------------------------- */
273 #ifdef __TI_COMPILER_VERSION__ // for TI compiler only
274 printf("Memory: %d bytes\n", &kernel_size);
275 #endif
277 printf("Cycles: %lld (order=%d) %lld (order=%d)\n",form_cycle[0],CYCLE_FORMULA_ORDER_PT1,form_cycle[1],CYCLE_FORMULA_ORDER_PT2);
279 }
280 /* ======================================================================= */
281 /* End of file: DSPF_dp_lud_d.c */
282 /* ------------------------------------------------------------------------*/
283 /* Copyright (c) 2013 Texas Instruments, Incorporated. */
284 /* All Rights Reserved. */
285 /* ======================================================================= */