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>
  43
  44 /* ======================================================================= */
  45 /* Interface header files                                                  */
  46 /* ======================================================================= */
  47 #include <ti/dsplib/dsplib.h>
  48
  49 /* ======================================================================= */
  50 /* define kernel size constant                                             */
  51 /* ======================================================================= */
  52 #if defined(__TI_EABI__)
  53 #define kernel_size _kernel_size
  54 #endif
  55
  56 /* ======================================================================= */
  57 /* extern parameters                                                       */
  58 /* ======================================================================= */
  59 extern char kernel_size;
  60
  61 /* ======================================================================= */
  62 /* define size constants                                                   */
  63 /* ======================================================================= */
  64
  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
  71
  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];
  78
  79 #pragma DATA_SECTION(U,".matrix")
  80 #pragma DATA_ALIGN(U,8);
  81 double U[MAX_U_SIZE];
  82
  83 #pragma DATA_SECTION(V,".matrix")
  84 #pragma DATA_ALIGN(V,8);
  85 double V[MAX_V_SIZE];
  86
  87 #pragma DATA_SECTION(U1,".matrix")
  88 #pragma DATA_ALIGN(U1,8);
  89 double U1[MAX_U_SIZE];
  90
  91 #pragma DATA_SECTION(diag,".matrix")
  92 #pragma DATA_ALIGN(diag,8);
  93 double diag[MAX_NROWS];
  94
  95 #pragma DATA_SECTION(superdiag,".matrix")
  96 #pragma DATA_ALIGN(superdiag,8);
  97 double superdiag[MAX_NROWS];
  98
  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;
 108
 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         }
 187
 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     }
 195
 196
 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     }
 211
 212     printf("test=%d:  %dx%d matrix has rank=%d\n",test,Nrows,Ncols,rank);
 213
 214   } /* for (test=1;test<=10;test++) */
 215
 216 }
 217
 218 /* ======================================================================= */
 219 /*  End of file:  DSPF_dp_svd_example.c                                    */
 220 /* ------------------------------------------------------------------------*/
 221 /*            Copyright (c) 2013 Texas Instruments, Incorporated.          */
 222 /*                           All Rights Reserved.                          */
 223 /* ======================================================================= */
 224