mexsvec_71.c

这事很不错的半定规划matlab软件,this is the important software of SDP.

/***********************************************************************
* mexsvec.c : C mex file 
*
*   x = mexsvec(blk,A,isspx,type); 
* 
*   Input: A     = nxn matrix
*          isspx = 0, store x as a dense vector
*                  1, store x as a sparse vector.
*          type = 0, stack upper triangular part of A column-wise 
*                 1, stack lower triangular part of A row-wise 
*
* SDPT3: version 3.0
* Copyright (c) 1997 by
* K.C. Toh, M.J. Todd, R.H. Tutuncu
* Last Modified: 2 Feb 01   
***********************************************************************/

#include <math.h>
#include <mex.h>

/**********************************************************
* single block: stack upper triangular part of A column-wise 
* into a column vector
**********************************************************/
void svec1(int n, double r2, 
           double *A, int *irA, int *jcA, int isspA, 
           double *B, int *irB, int *jcB, int isspB) 

{  int idx, i, j, jn, k, kstart, kend, idxj; 
   
   if (!isspB & !isspA) { 
      idx = 0; 
      for (j=0; j<n; j++) { 
          jn = j*n; 
          for (i=0; i<j; i++) { 
              B[idx] = A[i+jn]*r2; 
              idx++; } 
          B[idx] = A[j+jn]; 
          idx++; 
      }
   } 
   else if (isspB & !isspA) { 
      idx = 0; 
      idxj = 0; 
      for (j=0; j<n; j++) { 
          jn = j*n; 
          idxj += j; 
          for (i=0; i<j; i++) { 
              irB[idx] = i+idxj;               
              B[idx]   = A[i+jn]*r2; 
              idx++; } 
          irB[idx] = j+idxj; 
          B[idx]   = A[j+jn]; 
          idx++; 
      }  
      jcB[1] = idx;  
   }
   else if (!isspB & isspA) {
      idx = 0; 
      for (j=0; j<n; j++) { 
          idx += j; 
          kstart = jcA[j]; kend = jcA[j+1]; 
          if (kstart < kend) { 
             for (k=kstart; k<kend; k++) { 
                i = irA[k]; 
                if (i >= j) { break; } 
                B[idx+i] = A[k]*r2; 
             }
             if (i == j) {  B[idx+i] = A[k]; }        
	  }
      }
   } 
   else if (isspB & isspA) {
      idx = 0; 
      idxj = 0; 
      for (j=0; j<n; j++) { 
          idxj += j; 
          kstart = jcA[j];  kend = jcA[j+1];
          if (kstart < kend) { 
             for (k=kstart; k<kend; k++) { 
                i = irA[k]; 
                if (i >= j) {  break;  } 
                irB[idx] = i+idxj;              
                B[idx]   = A[k]*r2; 
                idx++; } 
             if ( i == j) { 
                irB[idx] = j+idxj; 
                B[idx]   = A[k]; 
                idx++;  }
          }  
      }
      jcB[1] = idx;  
   }
return;
}
/**********************************************************
* multiple sub-blocks: stack upper triangular part of A 
* column-wise into a column vector
**********************************************************/
void svec2(int n, int numblk, int *cumblksize, int *blknnz, 
           double r2, 
           double *A, int *irA, int *jcA, int isspA, 
           double *B, int *irB, int *jcB, int isspB) 

{  int idx, i, j, jn, l, jstart, jend, istart;
   int rowidx, idxj, k, kstart, kend; 
   
   if (!isspB & !isspA) { 
      idx = 0; 
      jstart = 0; jend = 0; 
      for (l=0; l<numblk; l++) { 
  	  jend = cumblksize[l+1];   
          istart = jstart; 
          for (j=jstart; j<jend; j++) { 
              jn = j*n; 
              for (i=istart; i<j; i++) { 
                  B[idx] = A[i+jn]*r2; 
                  idx++; } 
              B[idx] = A[j+jn]; 
              idx++;    }
          jstart = jend; 
      }
   }
   else if (isspB & !isspA) { 
      idx = 0; 
      jstart = 0; jend = 0; 
      for (l=0; l<numblk; l++) { 
  	  jend = cumblksize[l+1];   
          istart = jstart;  
          idxj = 0; 
          for (j=jstart; j<jend; j++) { 
              jn = j*n; 
              idxj += j-jstart;
              rowidx = blknnz[l]-istart+idxj; 
              for (i=istart; i<j; i++) { 
                  irB[idx] = rowidx+i;               
                  B[idx]   = A[i+jn]*r2; 
                  idx++; } 
              irB[idx] = rowidx+j; 
              B[idx]   = A[j+jn]; 
              idx++; 
          }
          jstart = jend; 
      }  
      jcB[1] = idx;  
   } 
   else if (!isspB & isspA) { 
      jstart = 0; jend = 0; 
      for (l=0; l<numblk; l++) { 
  	  jend = cumblksize[l+1];  
          istart = jstart;
          idx = blknnz[l]; 
          for (j=jstart; j<jend; j++) { 
              idx += j-jstart; 
              kstart = jcA[j]; kend = jcA[j+1]; 
              if (kstart < kend) { 
                 for (k=kstart; k<kend; k++) { 
                     i = irA[k]; 
                     if (i >= j) { break; } 
                     B[idx+i-istart] = A[k]*r2; 
                  }
                  if (i == j) {  B[idx+i-istart] = A[k]; }        
	      }
          }
          jstart = jend; 
      }  
   }
   else if (isspB & isspA) {
      idx = 0; 
      jstart = 0; jend = 0; 
      for (l=0; l<numblk; l++) { 
  	  jend = cumblksize[l+1];  
          istart = jstart;
          idxj = 0; 
          for (j=jstart; j<jend; j++) { 
              idxj += j-jstart;
              rowidx = blknnz[l]-istart+idxj;
              kstart = jcA[j];  kend = jcA[j+1];
              if (kstart < kend) { 
                 for (k=kstart; k<kend; k++) { 
                     i = irA[k]; 
                     if (i >= j) {  break;  } 
                     irB[idx] = rowidx+i;              
                     B[idx]   = A[k]*r2; 
                     idx++; } 
                 if (i == j) { 
                    irB[idx] = rowidx+j; 
                    B[idx]   = A[k]; 
                    idx++;  }
              }  
          }
          jstart = jend; 
      }
      jcB[1] = idx;  
   } 
return;
}
/**********************************************************
* single block: stack upper lowere part of A row-wise 
* into a column vector
**********************************************************/
void svec3(int n, double r2, 
           double *A, int *irA, int *jcA, int isspA, 
           double *B, int *irB, int *jcB, int isspB) 

{  int idx, rowidx, i, j, jn, k, kstart, kend; 
   
   if (!isspB & !isspA) { 
      idx = 0; 
      for (i=0; i<n; i++) { 
          for (j=0; j<i; j++) { 
              B[idx] = A[i+j*n]*r2; 
              idx++; } 
          B[idx] = A[j+j*n]; 
          idx++; 
      }
   } 
   else if (isspB & !isspA) { 
      idx = 0; 
      rowidx = 0; 
      for (i=0; i<n; i++) { 
          rowidx += i; 
          for (j=0; j<i; j++) { 
              irB[idx] = j+rowidx;               
              B[idx]   = A[i+j*n]*r2; 
              idx++; } 
          irB[idx] = j+rowidx; 
          B[idx]   = A[j+j*n]; 
          idx++; 
      }  
      jcB[1] = idx;  
   }
   else if (!isspB & isspA) {
      for (j=0; j<n; j++) { 
          kstart = jcA[j]; kend = jcA[j+1]; 
          for (k=kstart; k<kend; k++) { 
               i = irA[k]; 
               if (j < i) { 
 		  idx = i*(i+1)/2; 
		  B[j+idx] = A[k]*r2; }
               else if (j==i) {
		  idx = i*(i+1)/2; 
		  B[j+idx] = A[k]; }        
	  }
      }
   } 
   else if (isspB & isspA) {
      idx = 0; 
      for (j=0; j<n; j++) { 
          kstart = jcA[j];  kend = jcA[j+1];
          for (k=kstart; k<kend; k++) { 
              i = irA[k]; 
              if (j < i) {                  
                 irB[idx] = j + i*(i+1)/2;              
                 B[idx] = A[k]*r2;   
                 idx++; }
              else if (j==i) {
                 irB[idx] = j + i*(i+1)/2; 
                 B[idx] = A[k]; 
                 idx++;  
              }
          } 
      }
      jcB[1] = idx;  
   }
return;
}
/**********************************************************
* multiple sub-blocks: stack lower triangular part of A 
* row-wise into a column vector
**********************************************************/
void svec4(int n, int numblk, int *cumblksize, int *blknnz, 
           double r2, 
           double *A, int *irA, int *jcA, int isspA, 
           double *B, int *irB, int *jcB, int isspB) 

{  int idx, i, i1, j, l, jstart, jend, istart;
   int rowidx, idx2, k, kstart, kend; 
   
   if (!isspB) { 
      for (l=0; l<numblk; l++) { 
	  jstart = cumblksize[l]; jend = cumblksize[l+1];  
          istart = jstart;
          idx = blknnz[l]; 
          for (j=jstart; j<jend; j++) { 
              kstart = jcA[j]; kend = jcA[j+1]; 
              idx2 = idx + j-jstart; 
              for (k=kstart; k<kend; k++) { 
                  i = irA[k]; 
                  if (j < i) {
		     i1 = i-istart; 
                     B[idx2+i1*(i1+1)/2] = A[k]*r2; }
                  else if (j==i) {
		     i1 = i-istart; 
                     B[idx2+i1*(i1+1)/2] = A[k]; }        
              }
          }
      }  
   }
   else {
      idx = 0; 
      for (l=0; l<numblk; l++) { 
	  jstart = cumblksize[l]; jend = cumblksize[l+1];  
          istart = jstart;
          for (j=jstart; j<jend; j++) {               
	      rowidx = blknnz[l] + (j-jstart); 
              kstart = jcA[j];  kend = jcA[j+1];
              for (k=kstart; k<kend; k++) { 
		   i = irA[k]; 
                   if (j < i) {
                      i1 = i-istart; 
		      irB[idx] = rowidx + i1*(i1+1)/2;              
                      B[idx]   = A[k]*r2; 
                      idx++; } 
                   else if (j==i) { 
                      i1 = i-istart; 
                      irB[idx] = rowidx + i1*(i1+1)/2; 
                      B[idx]   = A[k]; 
                      idx++; }
	      }
          }
      }
      jcB[1] = idx;  
   }
return;
}
/**********************************************************
* Complex case: 
* single block: stack upper triangular part of A column-wise 
* into a column vector 
**********************************************************/
void svec1cmp(int n, double r2, 
           double *A, int *irA, int *jcA, int isspA, 
           double *B, int *irB, int *jcB, int isspB,
           double *AI, double *BI) 

{  int idx, i, j, jn, k, kstart, kend, idxj; 
   
   if (!isspB & !isspA) { 
      idx = 0; 
      for (j=0; j<n; j++) { 
          jn = j*n; 
          for (i=0; i<j; i++) { 
              B[idx]  = A[i+jn]*r2; 
              BI[idx] = AI[i+jn]*r2; 
              idx++; } 
          B[idx]  = A[j+jn]; 
          BI[idx] = AI[j+jn]; 
          idx++; 
      }
   } 
   else if (isspB & !isspA) { 
      idx = 0; 
      idxj = 0; 
      for (j=0; j<n; j++) { 
          jn = j*n; 
          idxj += j; 
          for (i=0; i<j; i++) { 
              irB[idx] = i+idxj;               
              B[idx]   = A[i+jn]*r2; 
              BI[idx]  = AI[i+jn]*r2; 
              idx++; } 
          irB[idx] = j+idxj; 
          B[idx]   = A[j+jn]; 
          BI[idx]  = AI[j+jn]; 
          idx++; 
      }  
      jcB[1] = idx;  
   }
   else if (!isspB & isspA) {
      idx = 0; 
      for (j=0; j<n; j++) { 
          idx += j; 
          kstart = jcA[j]; kend = jcA[j+1]; 
          if (kstart < kend) { 
             for (k=kstart; k<kend; k++) { 
                i = irA[k]; 
                if (i >= j) { break; } 
                B[idx+i]  = A[k]*r2; 
                BI[idx+i] = AI[k]*r2; 
             }
             if (i == j) {  B[idx+i] = A[k]; BI[idx+i] = AI[k];}        
	  }
      }
   } 
   else if (isspB & isspA) {
      idx = 0; 
      idxj = 0; 
      for (j=0; j<n; j++) { 
          idxj += j; 
          kstart = jcA[j];  kend = jcA[j+1];
          if (kstart < kend) { 
             for (k=kstart; k<kend; k++) { 
                i = irA[k]; 
                if (i >= j) {  break;  } 
                irB[idx] = i+idxj;              
                B[idx]   = A[k]*r2; 
                BI[idx]  = AI[k]*r2; 
                idx++; } 
             if (i == j) { 
                irB[idx] = j+idxj;