Fixed following bugs:
[ep-processor-libraries/dsplib.git] / ti / dsplib / examples / SVD_dp_rank_ex / DSPF_dp_svd_example.c
1 /* ======================================================================= */
2 /* DSPF_dp_svd_example.c -- Singular Value Decomposition example           */
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, INCSVDING, 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 (INCSVDING, 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 /*  (INCSVDING 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 <string.h>
42 #include <math.h>
44 /* ======================================================================= */
45 /* Interface header files                                                  */
46 /* ======================================================================= */
47 #include <ti/dsplib/dsplib.h>
49 /* ======================================================================= */
50 /* define kernel size constant                                             */
51 /* ======================================================================= */
52 #if defined(__TI_EABI__)
53 #define kernel_size _kernel_size
54 #endif
56 /* ======================================================================= */
57 /* extern parameters                                                       */
58 /* ======================================================================= */
59 extern char kernel_size;
61 /* ======================================================================= */
62 /* define size constants                                                   */
63 /* ======================================================================= */
65 #define MAX_NROWS   150 /* max is 150 */
66 #define MAX_NCOLS   150 /* max is 150 */
67 #define MAX_SIZE   (MAX_NROWS*MAX_NCOLS)
68 #define MAX_U_SIZE (MAX_NROWS*MAX_NROWS)
69 #define MAX_V_SIZE (MAX_NCOLS*MAX_NCOLS)
70 #define MAX_V_VECTOR_SIZE  MAX_NROWS
72 /* ======================================================================= */
73 /* double word alignment of arrays                                         */
74 /* ======================================================================= */
75 #pragma DATA_SECTION(A,".matrix")
76 #pragma DATA_ALIGN(A,8);
77 double A[MAX_SIZE];
79 #pragma DATA_SECTION(U,".matrix")
80 #pragma DATA_ALIGN(U,8);
81 double U[MAX_U_SIZE];
83 #pragma DATA_SECTION(V,".matrix")
84 #pragma DATA_ALIGN(V,8);
85 double V[MAX_V_SIZE];
87 #pragma DATA_SECTION(U1,".matrix")
88 #pragma DATA_ALIGN(U1,8);
89 double U1[MAX_U_SIZE];
91 #pragma DATA_SECTION(diag,".matrix")
92 #pragma DATA_ALIGN(diag,8);
93 double diag[MAX_NROWS];
95 #pragma DATA_SECTION(superdiag,".matrix")
96 #pragma DATA_ALIGN(superdiag,8);
97 double superdiag[MAX_NROWS];
99 /* ======================================================================= */
100 /* Main -- Top level driver for testing the algorithm                      */
101 /* ======================================================================= */
102 void main(void) {
103         
104   int test,Nrows,Ncols;
105   int row,col,status;
106   int rank,loop_limit;
107   double tolerance=0.000001;
109   /* --------------------------------------------------------------------- */
110   /* process test cases                                                    */
111   /* --------------------------------------------------------------------- */
112   for (test=1;test<=10;test++) {
113         switch(test) {
114         case  1: { Nrows=4;Ncols=2; // 4x2
115                        A[0]=1;A[1]=2;A[2]=3;A[3]=4;A[4]=5;A[5]=6;A[6]=7;A[7]=8;
116                    break;
117                  }
118         case  2: { Nrows=3;Ncols=3; // 3x3
119                A[0]=1;A[1]=1;A[2]=1;
120                A[3]=1;A[4]=2;A[5]=3;
121                A[6]=1;A[7]=3;A[8]=6;
122                break;
123                  }
124         case  3: { Nrows=3;Ncols=3; // 3x3
125                A[0]=1;A[1]=2;A[2]=3;
126                A[3]=4;A[4]=5;A[5]=6;
127                A[6]=7;A[7]=8;A[8]=0;
128                break;
129                  }
130         case  4: { Nrows=3;Ncols=3; // 3x3 singular
131                A[0]=0;A[1]=1;A[2]=2;
132                A[3]=3;A[4]=4;A[5]=5;
133                A[6]=6;A[7]=7;A[8]=8;
134                break;
135                  }
136         case  5: { Nrows=5;Ncols=4; // 5x4
137                A[ 0]=1;A[ 1]=0;A[ 2]=0;A[ 3]=0;
138                A[ 4]=0;A[ 5]=0;A[ 6]=0;A[ 7]=4;
139                A[ 8]=0;A[ 9]=3;A[10]=0;A[11]=0;
140                A[12]=0;A[13]=0;A[14]=0;A[15]=0;
141                A[16]=2;A[17]=0;A[18]=0;A[19]=0;
142                break;
143                  }
144         case  6: { Nrows=8;Ncols=5; // 8x5
145                A[ 0]=22;A[ 1]=10;A[ 2]= 2;A[ 3]=  3;A[ 4]= 7;
146                A[ 5]=14;A[ 6]= 7;A[ 7]=10;A[ 8]=  0;A[ 9]= 8;
147                A[10]=-1;A[11]=13;A[12]=-1;A[13]=-11;A[14]= 3;
148                A[15]=-3;A[16]=-2;A[17]=13;A[18]= -2;A[19]= 4;
149                A[20]= 9;A[21]= 8;A[22]= 1;A[23]= -2;A[24]= 4;
150                A[25]= 9;A[26]= 1;A[27]=-7;A[28]=  5;A[29]=-1;
151                A[30]= 2;A[31]=-6;A[32]= 6;A[33]=  5;A[34]= 1;
152                A[35]= 4;A[36]= 5;A[37]= 0;A[38]= -2;A[39]= 2;
153                    break;
154                  }
155         case  7: { Nrows=5;Ncols=5; // 5x5
156                A[ 0]=1;A[ 1]=0;A[ 2]=0;A[ 3]=0;A[ 4]=0;
157                A[ 5]=0;A[ 6]=1;A[ 7]=0;A[ 8]=0;A[ 9]=0;
158                A[10]=0;A[11]=0;A[12]=1;A[13]=0;A[14]=0;
159                A[15]=0;A[16]=0;A[17]=0;A[18]=1;A[19]=0;
160                A[20]=0;A[21]=0;A[22]=0;A[23]=0;A[24]=1;
161                    break;
162                  }
163         case  8: { Nrows=4;Ncols=5; // 4x5
164                A[ 0]=1;A[ 1]=0;A[ 2]=0;A[ 3]=0;A[ 4]=2;
165                A[ 5]=0;A[ 6]=0;A[ 7]=3;A[ 8]=0;A[ 9]=0;
166                A[10]=0;A[11]=0;A[12]=0;A[13]=0;A[14]=0;
167                A[15]=0;A[16]=4;A[17]=0;A[18]=0;A[19]=0;
168                    break;
169                  }
170         case  9: { Nrows=2;Ncols=4; // 2x4
171                    A[0]=1;A[1]=2;A[2]=3;A[3]=4;
172                    A[4]=5;A[5]=6;A[6]=7;A[7]=8;
173                    break;
174                  }
175         case 10: {
176                   Nrows=16;
177                       Ncols=Nrows;
178                       srand(1);
179                       for (row=0;row<Nrows;row++) {
180                         for (col=0;col<Ncols;col++) {
181                           A[row*Ncols+col]=(double)(rand())/((double)RAND_MAX);
182                         }
183                       }
184                   break;
185                 }
186         }
188     /* ------------------------------------------------------------------- */
189     /* decompose A into Q and R where A=U*D*transpose(V) using optimized C */
190     /* ------------------------------------------------------------------- */
191     status=DSPF_dp_svd(Nrows,Ncols,A,U,V,U1,diag,superdiag);
192     if (status==-1) {
193       printf("optimized C SV decomposition failed!\n");
194     }
197     /* ------------------------------------------------------------------- */
198     /* determine matrix rank by counting non singular values               */
199     /* ------------------------------------------------------------------- */
200     if (Nrows>=Ncols) {
201       loop_limit=Ncols;
202     } else {
203       loop_limit=Nrows;
204     }
205     rank=0;
206     for (row=0;row<loop_limit;row++) {
207        if (fabs(diag[row])>tolerance) {
208          rank++;
209        }
210     }
212     printf("test=%d:  %dx%d matrix has rank=%d\n",test,Nrows,Ncols,rank);
214   } /* for (test=1;test<=10;test++) */
218 /* ======================================================================= */
219 /*  End of file:  DSPF_dp_svd_example.c                                    */
220 /* ------------------------------------------------------------------------*/
221 /*            Copyright (c) 2013 Texas Instruments, Incorporated.          */
222 /*                           All Rights Reserved.                          */
223 /* ======================================================================= */